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Abstract 

This  study  investigated  the  effects  of  the  interaction  between  the  viscous  boundary  layer 
and  the  shock  wave  produced  by  a  Mach  10  inviscid  optimized  waverider.  An  implicit. 
Roe  flux-splitting  algorithm,  developed  by  WL/FIMM,  was  used  to  solve  the  flow  field. 
A  validation  for  the  inviscid  version  of  the  CFD  algorithm  was  accomplished  by 
comparing  the  numerical  data  produced  by  the  CFD  code  to  the  analytic  results  derived 
by  Rasmussen,  and  by  comparison  to  results  of  the  explicit  version  of  the  same  Roe 
flux-splitting  code.  The  computational  results  compared  favorably.  The  inviscid  case 
studied  using  the  implicit  code  produced  results  identical,  for  all  practical  purposes,  to 
those  of  the  explicit  code,  though  approximately  twice  as  quickly.  The  results  of  the 
viscous  flow  case  matched  well  with  the  results  predicted  by  theory.  The  lift  to  drag 
ratio  calculated,  5.74,  is  comparable  to  the  results  of  other  researchers. 


I.  INTRODUCTION 


1.1  Background 

The  waverider  concept  was  introduced  in  1958  by  Nonweiler  and  Hilton  (24), 
with  an  expanded  discussion  subsequently  published  by  Nonweiler  in  1959.  The  three- 
dimensional  body  Hilton  and  Nonweiler  examined  was  deemed  the  "caret  wing"  because 
of  the  distinctive  caret  shape  of  its  base  plane  (see  Figure  1-1). 


Figure  1-1,  Caret  Wing  Waverider 


A  waverider  is  a  lifting  body  which  is  derived  from  a  known,  analytic  flow  field, 
such  as  flow  over  a  two-dimensional  wedge  or  flow  around  a  slender  cone.  The  body 
is  designed  so  that  the  shock  is  attached  along  the  waverider’ s  entire  leading  edge,  thus 
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capturing  the  high  pressure  region  behind  the  shock  entirely  on  the  bottom  surface  of 
the  waverider.  The  top  surfoce  of  the  vehicle  is  typically  designed  to  be  completely 
parallel  to  the  freestream,  thus  creating  no  wave  drag  or  pressure  drag.  The  substantial 
difference  in  pressure  between  the  freestream  (upper)  and  compression  (lower)  surfaces 
of  the  waverider  generates  lift.  This  maximizes  the  lift  to  drag  (L/D)  ratio,  one  of  the 
primary  goals  of  high  speed  aircraft  designs. 

Hilton’s  and  Nonwieler’s  initial  waverider  work  was  quickly  followed  up  by 
others,  with  the  efforts  concentrating  on  two-dimensional  generating  flow  fields,  such 
as  flow  over  wedges.  Rasmussen  (26)  extended  the  waverider  concept  to  axisymmetric 
flow  fields  generated  by  slender  conical  bodies  using  Hypersonic  Small  Disturbance 
Theory  (HSDT)  to  analyze  flow  around  the  generating  conical  bodies  and  to  develop 
simplified  equations  which  describe  the  shape  of  the  waverider  based  upon  the 
generating  flow  field.  Rasmussen’s  work  in  the  1980’s  concentrated  primarily  on 
optimizing  waveriders  for  specific  flight  conditions,  with  Mach  number  and  generating 
cone  angle  being  the  dominant  parameters. 

Thus  far,  all  the  work  mentioned  has  been  inviscid,  i.e. ,  the  effects  of  viscosity 
were  neglected.  Waveriders  tend  to  suffer  from  viscous  effects,  as  they  have  a 
comparatively  low  volume  and  a  very  large  wetted  surface  area;  thus  much  research  has 
been  focused  on  examining  viscous  flow  over  waveriders.  Bowcutt,  et  al.  (10,  and  11), 
examined  viscous  optimized  waveriders,  optimizing  for  maximum  L/D,  by  means  of 
integral  boundary  layer  methods.  Corda  (12)  also  investigated  viscous  optimized 
waveriders  derived  from  a  flow  field  generated  by  a  power-law  minimum  drag  body. 
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where  a  reference  temperature  method  was  used  to  calculate  the  skin  friction.  At  high 
Mach  numbers,  Corda’s  method  gives  higher  L/D  values  than  Bowcutt’s,  although 
Bowcutt’s  method  does  better  than  Corda’s  at  Mach  numbers  less  than  10.  McLaughlin 
(21)  has  investigated  the  use  of  chemically  reacting  conical  flow  for  the  generating  flow 
Beld.  Vanmol  (33)  optimized  his  waverider  shape  for  maximum  L/D  using 
aerodynamic  heating  as  a  constraint.  Chang  (6)  investigated  the  effect  of  viscous 
boundary  layer  /  shock  interaction  on  waverider  optimization  through  the  use  of  the 
hypersonic  interaction  parameter  x  (defined  in  Chapter  5).  Takashima  (31)  validated 
the  integral  boundary  layer  method,  using  a  full  Navier-Stokes  solver,  by  investigating 
a  viscous  optimized  waverider  at  Mach  6. 

1.2  Objectives 

There  are  several  reasons  for  performing  this  thesis  research.  The  first  of  these 
reasons  is  that  the  determination  of  viscous  effects  on  waveriders  has  by  no  means  been 
fully  examined.  For  instance,  the  sharp  leading  edges  and  nose  common  to  inviscid 
waveriders  do  not  take  into  account  the  extreme  aerodynamic  heating  such  vehicles  will 
encounter.  While  Takashima  (31)  investigated  rounded  leading  edges,  it  was  only  for 
one  variety  of  waverider  over  a  narrow  range  of  flight  conditions.  This  thesis  will 
significantly  contribute  to  the  database  of  computational  solutions  to  waverider  designs 
by  examining  conically  derived  waveriders  with  both  sharp  and  blunted  leading  edges, 
including  a  detailed  examination  of  viscous  interaction  effects,  allowing  a  direct 
comparison  between  the  inversely  designed,  analytic,  inviscid  flow  field,  and  the 
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numerically  computed  inviscid  and  viscous  flow  fields.  Another  primary  reason  for  this 
thesis  is  to  provide  validation  results  for  the  Wright  Laboratory’s  high  speed,  implicit, 
finite-volume  Navier-Stokes  code  for  a  vehicle  of  comparable  shape  to  the  forebodies 
that  have  been  considered  for  the  National  Aerospace  Plane  (NASP).  This  code  has, 
to  date,  been  run  on  only  one  three-dimensional  configuration,  the  X-24C.  Specific 
objectives  this  thesis  will  accomplish  are: 

(1)  Develop  the  body  coordinates  of  an  inviscidly  derived  hypersonic  waverider, 
and  modify  it  to  round  the  leading  edges.  The  edges  will  be  rounded  by  separating  the 
freestream  surface  from  the  compression  surface  by  an  appropriate  amount,  and  then 
fitting  a  curve  between  the  edge  points  to  provide  a  smooth  leading  edge. 

(2)  Generate  a  three-dimensional  grid  suitable  for  capturing  the  expected  shock 
structure  and  the  boundary  layer  on  the  body. 

(3)  Apply  both  the  Euler  and  the  full  Navier-Stokes  versions  of  the  Wright 
Laboratory’s  three-dimensional,  implicit.  Roe  flux  difference  splitting  flow  solver 
developed  by  Gaitonde  (15)  to  the  hypersonic  waverider  and  grid  developed  in 
paragraphs  (1)  and  (2)  above. 

(4)  Compare  the  inviscid  computational  results  to  inviscid  analytical  results  (i.e. , 
Rasmussen  (25)  )  for  similar  configurations,  and  then  examine  the  effects  of  viscosity 
by  comparing  the  viscous  and  inviscid  comfRitational  solutions. 

1.3  Methodolt^ 

Rasmussen  (25)  and  (26),  details  the  waverider  design  methodology  applied 
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herein.  The  process  is  essentially  an  inverse  design,  where  the  waverider  shape  is 
derived  from  streamlines  in  a  previously  known  flowfreld.  The  vehicle  investigated  in 
this  effort  is  based  on  axisymmetric  hypersonic  flow  past  a  slender,  right  circular  cone. 
A  Cartesian  coordinate  system,  with  the  X-axis  along  the  cone  center,  the  Y-axis 
pointing  downward,  and  the  Z-axis  in  the  spanwise  direction  is  used.  Figure  1-2 
illustrates  the  Cartesian  coordinate  system  used  in  the  analysis. 


Since  Hypersonic  Small  Disturbance  Theory  (HSDT)  is  more  amenable  to 
spherical  coordinates,  the  Cartesian  system  described  above  was  transformed  to  a 
spherical  coordinate  system  with  the  waverider’s  nose  as  the  origin,  the  angle  B 
measured  from  the  X-axis,  and  <f),  the  azimuthal  angle  measured  from  the  Y-axis  in  the 
Y-Z  plane.  The  basic  parameters  that  determine  the  conical  flow  geometry  are  5,  the 
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cone  semi-angle,  jS,  the  shock  angle,  and  a  the  ratio  of  shock  angle  to  cone  semi¬ 
angle.  Figure  1-2  illustrates  the  spherical  coordinates  used  herein. 

As  mentioned  previously,  the  waverider  is  derived  using  HSDT.  The  process 
is  given  in  some  detail  in  chapter  2,  but  the  basic  method  used  to  derive  the  waverider 
is  as  follows: 

1 .  The  trailing  edge  of  the  freestream  surface  is  defined  with  a  four-term,  sixth- 
order  polynomial  given  by  Rasmussen  (26),  Eqn  (2-4).  This  effort  will  use  a  parabolic 
top  waverider,  which  uses  only  the  first  two  terms  of  the  polynomial. 

2.  The  flowfield  is  solved  analytically.  For  hypersonic  flow  around  a  slender 
cone,  an  analytic  solution  to  the  Taylor-MacColl  equation  (7),  is  obtained  using  the 
HSDT  technique. 

3.  Now  that  the  generating  flow  field  is  analytically  defined,  streamlines  parallel 
to  the  freestream  are  traced  upstream  from  the  freestream  trailing  edge  until  they 
intersect  the  shock  wave  produced  by  the  generating  cone,  as  shown  in  Figure  1-3. 
This  defines  the  entire  freestream  surfiice  and  the  leading  edge  of  the  waverider. 

4.  From  the  leading  edge  defined  in  step  3,  streamlines  are  traced  downstream, 
using  the  HSDT  relations,  until  they  cross  the  chosen  baseplane,  as  shown  in 
Figure  1-3.  These  streamlines  define  the  waverider’ s  compression  surfiice. 

Now  that  the  waverider  body  is  completely  defined,  a  computational  grid  is 
generated  to  define  the  domain  in  which  the  CFD  algorithm  will  be  applied.  The 
GRIDGEN  package  (30)  was  used  to  generate  the  grid.  Ideally,  the  grid  will  scale  with 
the  waverider  in  the  streamwise  direction,  so  that  the  shock  wave  will  be  captured  at 
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qjproximateiy  the  same  grid  location  for  each  cross  section.  The  scaling  also  preserves 
computational  resources,  which  would  otherwise  be  wasted  by  calculating  regions  of  the 
flow  which  should  remain  at  freestream  conditions.  However,  due  to  viscous  effects, 
it  may  not  be  possible  to  maintain  the  shock  at  the  same  grid  location,  as  viscous 
hypersonic  interaction  will  cause  the  shock  to  have  a  steeper  angle  than  predicted  by 
inviscid  flow  theory.  This  is  particularly  true  close  to  the  waverider’s  nose  where 
boundary  layer  growth  is  rapid,  and  the  shock/boundary  layer  interaction  is  strong. 

The  CFD  algorithm  used  to  obtain  a  computational  flow  solution  is  an  implicit, 
flux  splitting  Roe-averaged,  finite-volume  scheme.  The  Roe-averaging  scheme  is 
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particularly  desirable  in  this  case  due  to  its  shock  capturing  capability  and  stability  at 
high  Mach  numbers.  The  algorithm,  as  implemented  by  Gaitonde  (16),  is  also  quite 
robust  relative  to  grid  skewness  and  accurate  to  second  order  in  both  space  and  time. 
Numerical  values  of  the  velocity  components  (u,v,w\  Mach  number,  pressure  coefficient, 
lift,  wave  drag  and  viscous  drag  obtained  from  the  CFD  code  will  be  compared  to 
analytic  data  (25),  and  other  numerical  solutions. 
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IL  WAVERTOER  SURFACE  FORMULATION  AND  GRID  GENERATION 


2.1  Waverider  Surface  Formulation 

The  equations  defining  a  given  waverider  body  are  derived  in  detail  in  Appendix 
A.  Briefly,  the  waverider  is  described  by  the  following  equations,  with  the  coordinate 
systems  illustrated  by  Figure  1-2.  Except  where  noted,  the  following  equations  are 
derived  by  Rasmussen  (26).  To  begin,  the  freestream  surface  is  described  by: 

/e  =  rjim  (2-1) 

and  the  compression  surface  by; 

tiQ^  -  8^)’  =  r  (<i))(p'  -  by  (2-2) 

Equations  (2-1)  and  (2-2)  are  given  in  spherical  coordinates,  with  the  nose  of  the 
generating  cone  being  at  the  origin.  Note  that  r  =  r,((i>)  is  the  leading  edge  of  the 
waverider,  and  <|>  is  the  included  angle  measured  from  the  line  of  symmetry  to  the 
leading  edge.  Note  also  that  the  leading  edge  is  also  the  intersection  of  the  waverider 
body  with  the  shock. 

The  relation  of  the  cone  half  angle,  5,  to  the  shock  angle,  P,  is  provided  by 
HSDTas 


a  =  JL  = 


1*1  + 


(2-3) 


The  freestream  surface  of  the  waverider  is  defined  by  the  four  term,  sixth  order 
polynomial 
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(2-4) 


Y  =  R^  +  AZ^+BZ*  +  CZ^ 

where  A,  5,  and  C  are  constants  which  determine  the  curvature  of  the  freestream 
surface  of  the  waverider.  The  polynomial  given  by  Eqn  (2-4)  must  satisfy  two 
conditions 


Y  = 


:  Z  =  0 


(2-5) 


o  cos(t>,  :  Z  =  a  sin<|), 
where  the  angle  <{>/  is  the  included  angle  at  the  baseplane. 

At  the  baseplane,  the  compression  surface,  as  derived  by  Rasmussen  (26)  is 
deflned  by: 


RM  =  1  + 


0^-1 


V. 


(2-6) 


The  same  form  of  Eqn  (2-6)  holds  for  the  entire  length  of  the  waverider.  The 
included  angle,  <|>,  however,  decreases  with  decreasing  X.  The  slK>ck  relation  parameter, 
a,  scales  such  that  =  0  at  the  nose  coincides  with  Oo  =  /?„.  The  shock  attachment 
condition,  i.e.,  the  requirement  that  the  shock  must  attach  at  the  leading  edge  of  the 
waverider,  must  still  be  met  at  each  Y-Z  cross  section,  and  is  enforced  by 


Y  =  a^s<|)^  Z  =  a^in«|>^ 


(2-7) 


An  equation  for  results,  since  it  is  known  that  a  ranges  from  =  R„  to  o,  =  a. 
Thus 
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\^ere  XJl^  is  the  non-dimensionalized  length  of  the  waverider. 


Finally,  from  Stecklein  (28),  at  an  arbitrary  Y-Z  cross  section  is  given  by 


/ 

(  \ 

2 

f  \ 

X 

w 

R 

1  -  — 

+ 

0^-1 

1 

\  ' 

1  ^ ) 

“J 

1  J 

(2-9) 


However,  the  waverider,  as  developed  by  Rasmussen  (2S)  is  unsuited  to  viscous, 
hypersonic  flow  due  to  its  very  sharp  leading  edges.  Such  sharp  edges  magnify  already 
large  heat  transfer  rates.  Stecklein  (28)  developed  a  Fortran  code  to  generate  a 
waverider  body  using  the  equations  developed  by  Rasmussen  (25).  In  the  present  effort, 
Stecklein’s  Fortran  code  was  modified  to  generate  a  waverider  suitable  for  viscous, 
hypersonic  flow  as  follows: 

1.  The  waverider  is  translated  so  that  the  nose  of  the  waverider  is  at  the  origin. 

2.  The  compression  surface  was  then  moved  1  cm  in  the  positive  Y  direction 
to  provide  the  necessary  thickness  at  the  leading  edges  to  blunt  them. 

3.  Using  linear  interpolation,  the  X  coordinate  at  which  the  waverider  was  1  cm 
thick  (in  Z)  was  determined,  as  was  the  Y  coordinate  of  the  centerline  of  the 
compression  surface.  The  freestream  surface  (at  centerline),  of  course,  does 
not  vary  with  X. 

4.  Splines  were  fit  from  the  interpolated  point  to  a  boundary  point  imposed  .33 
cm  upstream  of  the  interpolated  point,  in  both  the  X-Y  and  X-Z  planes.  The 


.33  cm  value  was  determined  through  experimenting  with  various  values  for 
the  boundary  point.  The  ends  of  the  splines  were  clamped,  i.e.,  they  were 
required  to  match  the  slope  of  the  known  portion  of  the  waverider  and,  at  the 
boundary  point,  an  (arbitrary)  large  value. 

5.  Several  Y-Z  planes  were  determined,  using  the  splines  generated  in  step  4  as 
outer  bounds,  by  interpolation  of  data  from  the  first  known  plane  of  the 
waverider. 

6.  The  main  portion  of  the  waverider  had  splines  fit  in  the  Y-Z  plane  at  the 
leading  edges.  Again,  the  splines  were  clamped  and  an  arbitrary  boundary 
point  was  imposed.  The  splines  were  generated  first  from  the  freestream 
surface  to  the  boundary  point,  and  then  from  the  boundary  point  around  to  the 
compression  surface. 

In  this  effort,  28  points  were  added  at  each  plane  using  the  spline.  This  number 
of  points  gave  adequate  definition  for  the  model  without  making  the  model  database  file 
overly  large.  The  boundary  point  was  not  actually  added  to  the  waverider  model  as  it 
would  have  caused  a  sharp  point,  the  very  condition  the  spline  is  intended  to  eliminate. 

The  spline  is  given  by  a  cubic  polynomial; 

PH.)  =  (2-10) 

2  D 

where  the  are  constant  for  each  interval  /.  The  coefficients  can  be  solved  for  fairly 
easily  by  using  a  divided  difference  table  for  (see  Appendix  B). 

dc  Boor  (14)  gives  the  coefficients  of  Eqn  (2-10)  as 
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dPix)  , 

q^  =  _^=P/(x^=. 


'•'  2  At  *’'  ' 


(Ml) 


^  ^  P!"  (t,)  ^  ^,-«-^,.i-2[T,,T,.i]g 


4.1 


(At)^ 


Figure  2-1  and  Figure  2-2  show  the  geometry  obtained  using  the  cubic  splines 
at  the  nose,  and  Figure  2-3  shows  the  geometry  of  the  IS*  plane  from  the  nose.  The 
15*  plane  is  depicted,  as  it  is  the  first  plane  described  by  Stecklein’s  code,  with  the 
addition  of  the  cubic  spline  at  the  leading  edge.  Note  in  Figure  2-1  the  flat  spot  at  Y 
-  0.0  caused  by  leaving  the  boundary  point  out  of  the  model  definition  file. 


Figure  2-1,  X-Y  Plane  of  Waverider  Nose 
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Figure  2-2,  X-Z  Plane  of  Waverider  Nose 

Figure  2-4  is  a  three-dimensional  view  of  the  waverider  model,  showing  the 
results  of  the  whole  process  described  above.  The  modified  Waverider  code  is  in 
Appendix  C. 
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Figure  2-3,  15“*  Plane  from  Nose 
2.1.1  Program  WAVERIDER 

The  methods  described  by  Rasmussen  (25)  and  (26),  outlined  in  the  previous 
section  were  developed  into  a  Fortran  program  called  WAVERIDER  by  Stecklein  (28). 
Modifications  to  WAVERIDER  were  required  to  adapt  the  waverider  model  the  Fortran 
program  produced  to  one  suitable  for  viscous,  hypersonic  flow,  as  WAVERIDER  was 
originally  intended  to  produce  a  waverider  model  for  an  inviscid  flow.  The  modified 
code  defines  the  surface  of  a  parabolic-top,  inviscid  optimized  waverider,  just  as  written 
by  Stecklein  (28),  and  then  adapts  it  to  viscous  flow  by  blunting  the  leading  edges  and 
nose  using  the  methods  outlined  in  the  previous  section.  It  must  be  noted  that  the  code 
formulates  a  half-body  about  the  X-Y  plane  of  symmetry.  Building  the  full  body  would 
yield  profit  only  if  a  non-zero  yaw  was  considered  (in  which  case  it  would  be  required). 
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Figure  2-4,  Waverider  Model 

or  if  the  flow  were  to  become  asymmetric,  which  can  occur  for  low  speed  flow.  No 
such  cases  will  be  examined  in  this  effort. 

WAVERIDER  contains  the  parameters  which  govern  the  generation  of  the 
inviscid  optimized  waverider:  freestream  Mach  number,  generating  cone  half  angle  and 
length,  maximum  spanwise  sweep  angle,  and  the  type  of  parabolic  freestream  surface. 
The  design  parameters  used  here  are  given  in  Table  2.1.  The  surface  grid  dimensions 
are  given  by  the  integers  listed  in  Table  2.2.  Note  that  these  values  do  not  necessarily 
correspond  to  the  number  of  points  used  in  the  CFD  grid;  they  apply  only  to  building 
the  waverider  body  database  around  which  a  CFD  grid  is  later  constructed. 
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Table  2-1,  WAVERIDER  Parameters 


Parameter 

Value 

Meaning 

MACH 

10.0 

Freestream  Mach  number 

D 

5.5° 

Generating  cone  half  angle,  5 

L 

88.41675 

Generating  cone  length,  meters 

PHIL 

50.0° 

Maximum  spanwise  sweep  angle,  (|>b 

TYPE 

0 

0  indicates  parabolic  top  waverider 

Table  2-2,  WAVERIDER  Dimensions 


Parameter 

Value 

Meaning 

MCAP 

45 

Number  of  points  which  define  freestream  surface 
at  a  given  Y-Z  crossplane 

NCAP 

90 

Number  of  Y-Z  planes  defining  waverider 

INCR 

91 

Total  number  of  points  defining  waverider  surface 
at  a  given  Y-Z  crossplane 

2.1.2  Subroutine  WAVEBODY 

This  subroutine  takes  the  equations  derived  by  Rasmussen  (25),  outlined  above, 
calculates  the  waverider  body  shape  in  spherical  coordinates,  scales  the  equations,  and 
finally  transforms  the  spherical  coordinates  to  cartesian  coordinates.  The  primary  inputs 
to  this  subroutine  are  the  parameters  defined  in  the  main  program,  WAVERIDER.  The 
primary  outputs  of  this  program  are  NCAP  x  INCR  ordered  triples  (  (x,y,z)  coordinates) 
that  describe  the  waverider  in  cartesian  coordinates. 


2.1.3  Subroutine  OUTDATA 

The  subroutine  OUTDATA  has  three  primary  functions.  The  first,  as  suggested 
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by  the  name,  is  to  build  the  output  file  which  describes  the  waverider  surface.  The 
second  is  to  use  a  cubic  spline  and  linear  interpolation  techniques  to  build  a  blunt  nose. 
The  third  is  to  use  a  cubic  spline  fit  to  blunt  the  leading  edges. 

Primary  input  to  OUTDATA  are  the  coordinates  generated  by  WAVEBODY,  the 
absolute  value  of  the  slope  at  the  boundary  points,  and  the  number  of  points  to  use 
within  the  cubic  spline  fit.  The  latter  two  items  are  read  in  when  the  program  is 
executed. 


Table  2-3,  OUTDAT  Parameters 


Parameter 

Value 

Meaning 

DIM 

15 

Number  of  points  used  to  define  spline,  also  the 
number  of  planes  added  to  define  the  waverider’s  nose 

SLOPE 

m 

Absolute  value  of  the  slope  of  the  spline  at  the 
boundary  point 

2.1.4  Subroutine  GEOM 

This  subroutine,  provided  by  Beran  (8),  determines  the  spacing  of  both  the  Y-Z 
crossplanes  and  the  distribution  of  points  within  those  crossplanes.  Subroutine  GEOM, 
as  suggested  by  the  name,  determines  the  crossplane  spacing  via  a  geometric 
progression. 

2.1.5  Subroutine  CUBSPL 

CUBSPL,  developed  by  de  Boor  (14),  calculates  the  coefficients  in  Eqn  (2-10). 
Input  to  CUBSPL  consists  of  IBCl,  a  flag  that  sets  the  boundary  condition  used;  Xi,...,tj, 
the  Z-coordinates  of  the  endpoints  of  the  intervals;  gi(Tj), . ,g/Xj),  the  Y-coordinate  of 
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the  endpoints;  and  Si(Ti),...,Sj(Xj),  the  slopes  at  the  endpoints  of  the  intervals.  For  I 
intervals  there  will  be  I+l  of  each  of  the  latter  three  parameters. 

2.2  Grid  Generation 

The  GRIDGEN  software  package  (30)  was  used  to  generate  the  three-dimensional 
grid  around  the  waverider  body.  There  are  three  main  steps  in  generating  a  grid  with 
GRIDGEN; 

1.  GRIDBLOCK:  In  GRIDBLOCK,  the  boundaries  of  the  grid  are  defined,  as 
are  the  computational  coordinates,  i\,  and  and  their  associated  dimensions  using  a 
Silicon  Graphics  IRIS  4D  workstation.  For  the  waverider,  one  block  was  used,  although 
GRIDGEN  is  capable  of  handling  multi-block  configurations. 

2.  GRIDGEN2D:  In  GRIDGEN2D,  the  distribution  of  points  on  each  of  the  six 
faces  of  each  block  is  performed,  resulting  in  the  generation  of  surface  grids  on  each 
face  of  each  block.  This  step  is  also  performed  on  a  Silicon  Graphics  IRIS  4D 
workstation. 

3.  GRIDGEN3D:  GRIDGEN3D  is  the  final  step,  and  requires  a  CRAY  super¬ 
computer.  GRIDGEN3D  takes  the  results  of  GRIDGEN2D  and,  using  one  of  several 
solution  metlKids,  generates  the  volumetric,  three-dimensional  grid.  The  solution  is 
initialized  with  an  algebraic,  trans-finite  interpolation  method.  The  volume  grid,  once 
initialized,  may  be  further  refined  by  running  one  of  several  elliptic  solvers  for  a  given 
number  of  iterations. 

For  the  waverider  under  investigation  here,  grid  generation  was  performed  using 


19 


the  following  basic  procedure.  In  GRBDBLOCK,  as  mentioned  above,  a  single  block 
defining  the  grid  boundaries  and  computational  coordinates  was  generated.  The 
upstream  boundary  of  the  block  stands  off  .05  meters  from  the  nose  of  the  waverider 
(see  Figure  2-7)  ,  with  the  domain  being  an  elliptical  cone  S  m  in  diameter  at  the  nose 
and  25  m  in  diameter  at  the  waverider’s  baseplane,  in  the  Y  direction.  In  the  Z 
direction  the  elliptical  cone  was  5  m  in  radius  at  the  nose  and  12  m  in  radius  at  the 
baseplane.  Figure  2-5  and  Figure  2-9  illustrate  the  computational  domain  at  the 
baseplane  and  inflow  boundaries  respectively.  The  dimensions  were  chosen  to  make  the 
computational  domain  as  small  as  possible  without  having  the  imposition  of  the 
freestream  boundary  conditions  interfere  with  the  flow  solution.  The  ^  coordinate  points 
downstream  from  the  nose,  t]  is  normal  to  the  waverider  body,  and  ^  forms  a  right 
handed  system,  running  spanwise.  Grid  dimensions  were  set  to  61  x  91  x  101  in  the  4, 
T|,  and  ^  directions,  respectively. 

In  GRIDGEN2D,  grids  were  generated  on  each  face  of  the  block  created  in 
GRIDBLOCK.  The  most  interesting  was  face  5,  the  t)^„  face.  Face  5  was  first  split 
into  two  subfaces,  one  containing  the  6  points  between  the  nose  of  the  waverider  and 
the  upstream  block  boundary,  and  the  other  being  the  waverider  surface  itself  The 
waverider  surface  grid  was  generated  using  the  lines  on  database  networks  selection 
from  the  GRIDGEN2D  subface  shq)e  menu.  In  other  words,  the  surface  grid  was 
forced  to  conform  to  the  waverider  surface.  Choice  10  from  the  subface  interior  point 
initialization  (algebraic  solver)  menu,  interp  u.v  and fit  to  parametric  surface,  was  used 
to  generate  the  surface  grid. 
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Face  4,  the  face  or  waverider  baseplane  face,  required  some  special  attention 
as  well.  One  break  point  was  set  at  both  the  beginning  and  end  points  of  the  cubic 
spline  on  the  waverider  siuface.  This  put  49  grid  points  on  the  compression  surface,  12 
were  specified  on  the  cubic  spline,  leaving  40  on  the  freestream  surface.  Points  were 
clustered  to  the  leading  edge  at  the  waverider  surface,  and  away  from  the  centerline. 
On  the  outer  boundary  a  breakpoint  was  set  almost  directly  above,  though  a  little  out, 
from  the  leading  edge.  Again,  49  points  were  specified  for  the  compression  region. 
Points  on  the  outer  boundary  were  clustered  away  from  the  breakpoint,  and  from  the 
centerline,  leading  to  a  clustering  near  the  middle  of  the  two  subedges.  This  spacing 
accomplished  two  things.  First,  it  kept  gridlines  from  overlapping  each  other,  or  cutting 
through  the  waverider  body,  and  secondly,  it  provided  for  grid  clustering  in  the  region 
where  the  shock  is  expected  to  form.  Figure  2-5  and  Figure  2-6,  on  the  following 
pages,  illustrate  the  results  of  this  process.  Figure  2-6  illustrates  the  grid  geometry  close 
to  the  surface  of  the  waverider. 
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Figure  2-6,  Baseplane  Grid,  Zoom  of  Cusp 


Similar  measures  were  taken  at  the  nose,  or  freestream  plane,  which  is  face  3, 
the  face.  Choice  1  from  the  subface  interior  point  initialization  (algebraic  solver) 
menu,  arclength  based  (Soni)  TFI  was  used  to  initialize  this  face,  and  all  other  faces 
except  face  5. 

On  faces  1  and  2,  the  and  faces,  GRIDGEN2D’s  elliptic  solver  was  used 
to  reduce  grid  skewness  near  the  waverider’s  nose.  Thomas-Middlecoff  control 
functions,  with  a  relaxation  factor  of  .25,  were  used  for  the  smoothing.  Thomas- 
Middlecoff  control  functions  were  used  as  they  provided  the  best  results.  The  small 
relaxation  factor  was  chosen  so  that  the  changes  to  the  grid  would  be  slow  enough  that 
when  the  grid  was  acceptably  smooth  the  process  could  be  stopped.  Approximately  200 
iterations  were  required  to  produce  an  acceptable  grid. 

The  rest  of  the  faces  of  the  block  were  rather  straightforward,  requiring  little 
special  treatment.  Points  were  clustered  toward  the  leading  edges,  the  nose,  and  the 
waverider  surface  to  improve  flow  resolution,  particularly  in  the  boundary  layer  region. 
Normal  to  the  surface,  the  first  point  is  spaced  .5  millimeter  above  the  surface.  This 
spacing  was  arrived  at  by  looking  at  a  converged  solution  from  a  much  coarser  grid. 
The  coarse  grid’s  first  point  was  .5  centimeters  above  the  surface.  At  the  baseplane 
there  were  about  three  points  in  the  boundary  layer,  so  the  spacing  in  the  fine  grid  was 
made  one  tenth  of  that,  to  provide  approximately  20  points  in  the  boundary  layer  at  the 
baseplane.  The  following  figures  illustrate  the  grids  produced  by  GRIIXjEN2D  on  the 
various  faces  of  the  block. 
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Figure  2-7,  Centerline  Grid,  Faces  1  &  2 


In  GRIDGEN3D,  as  mentioned  previously,  the  volume  grid  was  initialized  using 
an  algebraic  trans-flnite  interpolation  method.  This  produced  385  skewed  volumes  and 
19  negative  volumes.  GRIDGEN3D’s  elliptic  solver,  using  Laplace  control  functions, 
was  then  run  for  12  iterations  to  smooth  the  grid.  The  negative  volumes  were 
eliminated,  and  the  number  of  skewed  volumes  remained  the  same.  GRIDGEN3D 
required  approximately  5  minutes  to  run  on  the  CRAY  Y-MP  from  the  time  the  batch 
job  was  submitted  until  it  was  completed.  The  figures  depicting  the  grid  at  the  various 
boundaries  shown  on  the  previous  pages  are  the  output  of  GRIDGEN3D,  with  ghost 
points  (points  within  the  body  surface)  required  by  the  computational  solver  added. 

The  ghost  points  mentioned  above  are  required  because  finite  volume  codes,  such 
as  the  one  used  here,  solve  for  flux  across  cell  faces,  which  requires  knowing  the  flow 
at  cell  centers,  rather  than  at  the  grid  points.  Ghost  points  are  added  to  all  boundaries, 
so  that  cell  centers  can  be  calculated.  Figure  2-1 1  illustrates  the  addition  of  ghost  points 
to  a  simple  grid,  and  the  resulting  network  of  cell  centers. 
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Figure  2-11,  Addition  of  Ghost  Points  to  Grid 
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in.  NAVBER-STOKES  IMPLiaT  FLUX  SPLITTING  ALGORITHM 


3.1  Governing  Equations 

The  Navier-Stokes  equations  are  well  known,  and  many  sources  provide 
derivations  of  them.  Anderson  (3),  among  others,  lists  the  Navier-Stokes  equations  as 


Conllnuiv  SP  ,  ^  ^  ^  =  0  (3-1) 

dt  dx  dy  dz 

X-Momeruum  £21  =  -±  *  ^  ^  ^  (3-2) 

Dt  dx  dx  dy  dz 

Y-Momenlum  =  -  JP  .  (3-3) 

Dt  dy  dx  dy  dz 


Z-Momentum 


Dpw  ^  _dp  ^  ^ 

Dt  dz  dx  dy  dz 


dx.. 


dx. 


(3-4) 


j.  DpE,  d 

Energy  __  =  +  _ 

Dt  dx 


[kE] 

.1. 

\kE\ .  ± 

fk^] 

dx) 

dy 

[  dy)  dz 

1 

dx  dy  dz  dx  dy 

dz  dx  ^  dz 


(3-5) 
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where  E,  is  total  energy. 

The  shear  terms  are  defined  as: 


X,.  =  5,/(V-F)  + 


du  du 


dx  dx 


J  J 


(3-6) 


where  5  is  the  Kronecker  delta,  and  the  subscripts  correspond  to  standard  indicial 
notation.  Note  that  the  shear  stress  terms  are  symmetric,  i.e.  x^  = 

Anderson,  et  al.  (1)  presents  equations  (3-1)  through  (3-5)  in  vector  form  as 


where 


dU  dE  ^  dF  ^  dG 
dt  dx  dy  dz 


(3-7) 


U  = 


P 

P« 

pv 


pw 


(3-8) 


E  = 


P« 

P«^  +/»  -  \ 


pttv  -  X 


*y 


puw  -  X^ 

(£,  +  p)u  -  «x„  -  vx^  -  wx„  + 


(3-9) 
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F  = 


pv 

puv  +  p  - 


pv'  -  t 


yy 


pvw  -  X 


yi 


(£,  +  p)v  -  ux  -  vx  -  wx 

\\  t  t'y  xy  yy  yi 


G  = 


pw 

puw  +  p  -  x^ 


pvw  -  X 


>« 


pw^  - 
[(£,  +  p)w  -  ux^  -  VT^ 


+ 


(3-10) 


(3-11) 


The  form  of  the  Navier-Stokes  equations  given  in  equations  (3-7)  to  (3-11)  is 
much  more  convenient  than  that  in  equations  (3-1)  through  (3-5)  for  purposes  of  coding 
the  Navier-Stokes  equations  into  a  CFD  algorithm.  The  form  of  equation  (3-7)  is  easier 
to  code  than  that  in  equations  (3-1)  through  (3-S)  because  it  is  more  concise. 

Transforming  the  Navier-Stokes  equations  to  the  computational  coordinate  system 
is  necessary  so  that  the  equations  are  consistent  with  the  computational  space.  The 
general  transformed  form  of  the  Navier-Stokes  equations  is  (1): 
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(3-12) 


*  *  ’’■p'  ■  ^->1} 

*  -  ^-) "  -  ■^■) "  -  ^J]} '  “ 

The  Jacobian,  J,  is  the  matrix  of  inverse  metrics: 

j  ^  0  (3-13) 

3(*,y,^) 

See  Beran  (8)  for  the  complete  statement  of  the  Jacobian  matrix. 

Note  in  equation  (3-12)  that  the  flux  vectors  E,  F  and  G  have  been  separated  into 
inviscid  and  viscous  terms.  The  flux  vector  for  instance,  is  separated  into  Ei  and  E„ 
where: 

pu 

pu^  +  p 

E.  =  puv  £  = 

puw 

(E,  +  p)u 

The  viscous  terms  in  F  and  G  are  separated  from  the  inviscid  terms  in  a  similar  fashion. 

3.2  Discretization 

The  code  developed  at  WL/FIMM  and  ^plied  in  the  present  research  uses 
several  methods  in  conjunction  with  each  other  to  discretize  the  Navier-Stokes  equations. 
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The  viscous  terms  in  the  flux  vectors  are  handled  by  simple  centered  differences.  The 
inviscid  terms  are  handled  by  a  combination  of  flux  vector  splitting  and  flux  difference 
splitting  as  described  below. 

Consider  the  one-dimensional,  inviscid  model  equation; 


dt  dx 


(3-15) 


First  order  implicit  discretization  of  equation  (3-15)  with  the  backward  Euler 
method  produces 


ur  -  u” 

A? 


■*"^1 4  *  ■*^1  •  * 

- Z _ L  =  0 

Ax 


(3-16) 


where  the  subscript  i+1/2  denotes  the  interface  between  nodes  /  and  i+1. 
Linearization  of  a  typical  flux  term  yields 


where 


£"  V  =  £ " .  + 

'  T  '*T 

=  E” .  + 


At 


(  dE  dU') 


\^dU 


At 


(3-17) 
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(3-18) 


A  = 


dU 


8C/"  =  U”*^  -  U" 


Substituting  equation  (3-17)  into  equation  (3-16)  yields: 


_ 1  +  _ I _ 1 

At  Ax 


(3-19) 


Flux  vector  splitting,  in  the  form  of  the  Steger-Warming  algorithm  (29),  is  used 
on  the  left  hand  side  of  equation  (3-19),  while  Roe  flux  difference  splitting  (8),  in  the 
form  of  the  MUSCL  scheme  (37),  is  used  on  the  right  hand  side  of  equation  (3-19). 

The  three-dimensional  system  in  equation  (3-7)  is  discretized  very  much  like  the 
one-dimensional  model  system  above.  In  order  to  compute  the  right  hand  side,  or 
residual,  of  the  three-dimensional  version  of  equation  (3-19),  the  fluxes  in  each  direction 
are  successively  balanced.  Thus 


LHS^  =  _Ll_ 
LHS^  =  LHS^  + 

LHS,  =  LHS,  + 


Ax 

ff, .,(£/■) 

T _  T 

I  J 

f G,,.(C/")  -  G,. .(t/-)' 


(3-20) 


where  LHS|  refers  to  the  left  hand  side  of  the  three-dimensional  version  of  equation 
(3-19). 
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The  left-hand  side  of  equation  (3-20)  is  calculated  using  Steger-Warming  flux- 
vector  splitting  (29),  discussed  in  the  following  section.  Roe’s  upwind  method  (8)  is 
used  on  the  right  hand  side  of  equation  (3-20)  to  calculate  and  For 

example; 


^'•4  =  1 


E(U,)  *  E(U,)  -  -  t4) 


(3-21) 


The  terms  Ur  and  C4  obtained  from  the  MUSCL  approach  described  in  Section  3.4. 
The  indicates  a  Roe  averaged  term,  also  described  in  Section  3.4 


3.3  Flux  Vector  Splitting 

Flux  vector  splitting  is  a  well  known  technique  used  to  improve  the 
computational  efficiency  of  finite-difference  schemes,  as  well  as  to  make  the  schemes 
somewhat  more  robust.  The  basic  aim  is  to  split  the  flux  vectors  such  that  an  upwind 
finite-difference  scheme  may  be  used  at  all  points  within  the  flow.  This  is  done  by  the 
simple  expedient  of  separating  the  flux  vectors  into  upwind  and  downwind  parts.  In 
WL/FIMM’s  high  speed,  implicit,  flux-difference  splitting  CFD  code  is  used  on  the  left- 
hand  side  of  the  finite  difference  equation,  i.e.  the  left  hand  side  of  equation  (3-20). 
Steger  &  Warming  (29)  introduce  flux  vector  splitting  in  some  detail.  Briefly,  for  the 
three-dimensional  Navier-Stokes  equations,  flux-vector  splitting  proceeds  as  follows. 

Equation  (3-7),  neglecting  viscous  terms,  can  be  rewritten  as 


^  +  A— 
dt  dx 


dy  dz 


(3-22) 
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where  A,  B,  and  C  are  the  Jacobian  matrices 


dK 

—L  B  =  — i 
dU  dU 


(3-23) 


which  correspond  to  the  inviscid  portions  of  the  flux  vectors,  Fp  and  G,.  Note  that 
Ei=AU,  Fi=BU,  and  Gi=CU  by  the  first  order  homogenous  property  of  the  Euler 
equations. 

A  finite-difference  form  of  equation  (3-22)  is: 


w;  ^  d[a;w”]  ^  d[b;w;)  ^  d[c”w;) 

Ai  Ax  Ay  Az 

DE”  df;  dg; 

=  -  _  +  _  +  _ 

_  Ajc  Ay  ^  _ 

where  5  and  D  are  the  time  and  space  difference  operators,  respectively; 


(3-24) 


(3-25) 

(3-26) 

Now,  define  the  matrix  P  as  a  linear  combination  of  A,B,  and  C: 


P  =k^A  +  +  kf  (3-27) 

A^ere  the  k,  are  arbitrary  real  constants. 

The  system  in  equation  (3-22)  is  hyperbolic  if  there  exists  a  similarity  transformation 
such  that 
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Q-^PQ 


X,  0  0  0  0 
0  ^2  0  0  0 
0  0  Xj  0  0 
0  0  0  0 
0  0  0  0  x. 


=  A 


(3-28) 


The  (general)  flux  vector  is  given  by 


QKQ-^U  (^29) 

Note  that  the  Jacobian  can  be  diagonalized  due  to  the  hyperbolic  nature  of  the  inviscid 
fluxes.  The  eigenvalues,  X;,  of  the  5x5  matrix.  A,  are  arbitrary,  but  real.  The  matrix  Q 
and  its  inverse  are  calculated  as  follows: 


Q=MT 

The  matrices  M,  T,  and  their  inverses  are  given  by  Warming  (35). 

The  splitting  comes  in  through  the  eigenvalue  matrix,  A.  It  is  split  into  two 
matrices,  A^  and  A'  via  the  splitting  of  the  eigenvalues,  X;: 


V  =  *  i\i)  V  =  - 1^,1) 

Thus,  the  flux  vector  is  split  by  using  A*  and  A'  in  place  of  A  in  equation  (3-29), 


so: 


^  =  QA*Q  ‘f/  ^  =  QA  Q  -‘t/  (^32) 

Hence,  in  equation  (3-24)  the  terms  operated  on  by  the  D/Ax  operator  can  be  rewritten 


as: 
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(3-33) 


T 


=  a;sc/,  t  -  (a,:,su,.,  *  a;w} 


3.4  Roe  Algorithm 

Beran  (8),  presents  the  Roe  scheme  for  a  one  dimensional  system  of  equations 
such  as  in  equation  (3-1 S),  as 

=  u;  -  -  i/".) 

(3-34) 


As  mentioned  in  Section  3.2,  the  code  developed  by  WL/FIMM  uses  a  variation 
of  the  Roe  scheme.  Monotonic  Upstream  Scheme  for  Conservation  Laws,  or  MUSCL. 
MUSCL  is  a  form  of  a  TVD,  total  variation  decreasing,  scheme.  Yee  (37)  presents  the 
MUSCL  scheme  using  a  simple,  one  dimensional  hyperbolic  equation.  Starting  from 
equation  (3-34)  the  MUSCL  scheme  replaces  the  1/^.1  and  Ui  terms  with  and  lJ^,a 
respectively.  and  are  given  b> 


=  (/, 


/♦I 


(3-35) 


and 
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where: 


t/,‘,  =  U,  -  ^(l  -  ^A,.,  +  (l  +  ^)a,.. 


(3-36) 


11  = 


-1 :  upwind  scheme 

0:  Fromm  scheme 

1/3;  third  order  upwind  biased  scheme 

1;  three  point  central  difference  scheme 


A.,,(/  =  minmo<{u.^,  -  U.,(o{U.  -  (/..,)) 


(3-37) 


(3-38) 


and,  finally 


=  minmoc(u.^^  -  U.,ay(U.^^  -  C/,,))  (3-39) 


1  ^  ^  3  -  ii 

1  <  o  ^ - J. 

1  -  Ti 


ri  ^  1 


(3-40) 


minmod(x,y)  =  5gn(x)Tnax{o,min[|ac|,yjg77(ar)]}  (3-41) 

Note  that  sgn(x)  means  the  sign  of  the  variable  x.  Note  also  that  the  minmod  slope 
limiter,  equation  (3-41),  is  not  the  only  limiter  that  may  be  applied. 

The  in  the  above  equations  refers  to  Roe  averaging,  which  is  given  by: 


u  = 


Pr  ^r 

(3-42) 


The  L  and  R  subr.cripted  variables  in  equation  (3-42)  refer  to  the  components  of  the  like 


41 


subscripted  vector  U  in  equations  (3-35)  and  (3-36). 


3.5  Boundary  and  Initial  Conditions 

There  are  four  basic  boundary  conditions  imposed  in  the  waverider  problem: 
freestream  conditions,  surface  conditions,  symmetry  conditions  and  no-change 
conditions.  The  initial  conditions  were  simply  uniform  flow  at  freestream  conditions. 

The  freestream  conditions  were  imposed  on  the  face  and  the  face. 
Fluxes  on  the  freestream  boundaries  were  calculated  using  the  values  for  Mach  number, 
temperature  and  pressure  assigned  in  the  initial  conditions.  On  the  face,  the  inflow 
condition  is  specified  by 

f/,_.  =  U,  (3-13) 

where  U  is  the  vector  given  in  equation  (3-8). 

Similarly,  the  freestream  condition  is  given  on  the  face  by 

(3-44) 

The  boundary  conditions  at  the  waverider’s  surface  are  a  bit  more  complicated. 
The  viscous  terms  are  calculated  directly  from  the  conditions  specified  at  the  surface  for 
temperature  and  zero  velocity.  Because  the  computational  grid  is  cell  centered,  the 
conditions  at  the  surface  are  imposed  by  setting  the  conditions  at  the  ghost  points  just 
within  the  surface.  Thus,  for  the  u-component  of  velocity 
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(3-45) 


identical  relations  are  used  for  the  other  velocity  components,  v,  and  w.  Furthermore, 
pressure  is  held  constant  between  the  ghost  point  and  the  first  point  above  the  surface; 


P  =  P 

i,l,*  ^  i,2,* 

Temperature  at  the  ghost  point  is  calculated  using; 


(3-46) 


7',.a  =  2  7-^,-  r  ^  (3-47) 

The  inviscid  fluxes  are  calculated  directly  from  the  pressure  and  the  metrics. 
Since  the  inviscid  energy  and  mass  fluxes  normal  to  the  surface  are  both  zero,  the  fluxes 


are  given  by; 


(3-48) 


Recall  from  above  that  the  pressure  is  assumed  constant  between  the  ghost  point  and  the 
first  point  above  the  surface.  Note  also  that  the  metrics  are  calculated  at  the  wall,  since 
wall  points  are  actually  grid  points. 

The  symmetry  boundary  condition  is  imposed  on  the  and  faces.  Flux 
in  the  ^  direction  is  set  to  zero  by  setting  the  flux  at  the  ghost  points  equal  to  the 
negative  of  the  flux  at  the  first  point  out  from  the  symmetry  plane.  For  the  face  this 
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is: 


P^v..  =  -P^v.2 

with  the  exact  same  form  for  the  face. 

A  no-change  condition  is  applied  to  the  outflow  boundary,  the  face.  The  no¬ 
change  condition  is  imposed  by  simply  setting  the  fluxes  at  the  ghost  point  equal  to 
those  at  the  first  point  upstream  of  the  boundary,  thus; 


(3-50) 


When  calculating  quantities  such  as  the  lift  to  drag  ratio,  L/D,  the  calculations  were 
simply  terminated  at  the  outflow  boundary. 


3.6  Computer  Code  Description 

Wright  Laboratory’s  three-dimensional,  implicit,  flux-splitting,  Navier-Stokes, 
Fortran  code  developed  by  Gaitonde  (15)  under  contract  to  WL/FIMM,  was  used  to 
solve  the  viscous,  hypersonic  flow  over  a  Mach  10,  inviscid-optimized  waverider, 
modified  for  viscous  flow.  The  code  was  debugged  on  the  ASD  CRAY  X-MP  super¬ 
computer  using  a  Silicon  Graphics  IRIS  workstation  as  an  interface.  The  debugged  code 
was  transferred  to  an  IRIS  workstation  in  the  AFIT  computer  laboratory  from  which  it 
was  sent  to  the  CRAY  Y-MP  at  the  Ohio  Supercomputer  Center  (OSC),  w^ere  all 
solution  runs  were  performed.  Access  to  the  OSC  CRAY  was  through  a  grant  for 
research  on  Numerical  Solution  oflnviscid/Viscous  Hypersonic  Flow  Around  a  Conically 
Derived  Waverider,  grant  number  PIS009-2.  Connection  to  the  OSC  CRAY  was  made 
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via  telnet  from  the  AFIT  IRIS  workstations. 


Two  data  files  were  required  for  execution  of  the  code:  cnldat  and  cnJgrd.  The 
first  file,  cnldat,  contains  information  such  as  the  version  of  the  code  to  run  (Navier- 
Stokes  or  Euler),  whether  to  use  an  implicit  or  an  explicit  formulation,  initial  conditions, 
boundary  conditions,  and  so  forth.  The  cnldat  file  has  essentially  four  parts:  solution 
integration  parameters,  flow  field  conditions,  boundary  conditions,  and  control 
parameters.  A  listing,  by  category,  of  the  input  data  file  follows.  A  sample  cnldat  file 
can  be  found  in  Appendix  D. 

The  first  table.  Table  3-1,  lists  the  parameters  that  control  the  implementation  of 
the  solution,  particularly  such  things  as  which  finite  difference  scheme  to  use,  i.e..  Roe, 
Lax-Wendroff,  or  Van  Leer.  Other  parameters  include  the  ending  iteration  number, 
whether  the  flow  model  is  Euler  or  Navier-Stokes,  what  the  CFL  stability  criteria  is,  and 
how  it  is  handled.  Table  3-2  details  the  flow  field  conditions,  particularly  the  ffeestream 
conditions  and  the  orientation  of  the  waverider,  along  with  some  data  about  the 
waverider  itself,  such  as  surface  temperature.  Table  3-3  lists  the  boundary  condition 
parameters,  and  the  node  range  in  which  they  apply. 


45 


Table  3-1,  Solution  Integration  Parameters 


ICASE 

Finite  Difference  Scheme:  l=Roe 

MEND 

Number  of  iterations 

INS 

Flow  model;  l=Navier-Stokes,  0=Euler 

IL,JL,KL 

Grid  Dimensions 

CFLMAX 

Maximum  CFL  number  allowed 

CFL 

Starting  CFL  number 

CFLEXP 

Number  of  iterations  before  CFL  doubles 

ICFL 

Number  of  iterations  between  CFL  increases 

IPC 

Driver:  l=Backward  Euler,  l=Predictor-Corrector 

IMPLT 

Implicit  vs.  Explicit:  0=Explicit,  l=Implicit 

Table  3-2,  Flow  Field  Conditions 


lADBWL 

Flag  to  indicate  adiabatic  wall  BC 

ALPHA 

Angle  of  attack 

PHI 

Yaw  angle 

TWALL 

Wall  temperature 

RMINF 

Freestream  Mach  number 

REL 

Reynolds  number  per  foot 

RLSCL 

Scaling  factor 

TINF 

Freestream  Temperature 

Table  3-3,  Boundary  Conditions 


me 

Start  and  end  points  for  a  ^  BC 

me 

Start  and  end  points  for  an  t|  BC 

KBC 

Start  and  end  points  for  a  ^  BC 

mCTYPE 

Boundary  condition  type  imposed 
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Table  3-4,  Control  Parameters 


IREAD 

0=Deadstart,  l=Restart 

IGRID 

Grid  format 

IP3DOP 

PlotSD  output  file  format 

MODPR 

Iterations  between  printing  convergence  data 

The  file  cnlgrd  contains,  in  binary  format,  the  grid  generated  using  GRIDGEN 
(as  described  in  Chapter  2).  The  GRIDGEN  output  was  converted  to  a  cell  centered 
grid  by  the  routine  REDUCE.F,  supplied  by  Gaitonde.  The  coordinate  axes  defined  in 
GRIDGEN,  ^  =  streamwise  axis,  r\  =  surface  normal  axis,  and  Q  =  spanwise  axis,  were 
maintained. 
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IV.  COMPUTATIONAL  RESULTS 


4.1  Inviscid  flow  Results 

The  design  point  for  the  inviscid  optimized  waverider  being  studied  is  Mach  10 
at  100,000  feet,  identical  to  that  studied  by  Stecklein  (28).  Table  4.1  lists  the  design 
parameters.  The  design  point  should,  for  the  inviscid  case,  produce  the  best  lift  to  drag 
(L/D)  ratio  possible  for  the  waverider  as  ftiat  is  the  condition  the  waverider  was 
optimized  for. 


Table  4-1,  Waverider  Design  Parameters 


Parameter 

Value 

Mach  Number 

10 

Altitude 

100,000  ft 

Freestream  Temperature 

406.7  °R 

Wall  Temperature 

530  °R 

Freestream  Density 

3.32  E-05  slugs/ft^ 

Freestream  Pressure 

22.43  Ibs/ft^ 

4.1.1  Unmodifled  Waverider 

In  order  to  establish  consistency  between  the  implicit  version  of  the  CFD 
algorithm  and  the  explicit  version,  an  inviscid  case,  using  a  waverider  model  and  grid 
identical  to  that  used  by  Stecklein  (28),  was  run  to  a  fully  converged  solution. 
Convergence  is  defined  as  occurring  when  the  residual  drops  below  10'^.  The  residual 
is  given  by 

For  plotting  purposes,  the  residual  as  defined  in  equation  (4-1)  was  normalized  by  the 
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Stecklein  (28)  computed  a  value  of  8.18812  for  L/D  in  his  work,  using  the 
explicit  version  of  the  code.  The  implicit  version  of  the  code  returned  a  value  of 
8.18813.  The  two  values  differed  in  the  fifth  decimal  place. 

The  main  difference  between  the  results  from  the  implicit  and  explicit  versions 
of  the  code  is  the  convergence  history.  The  implicit  case  was  fully  converged,  using 
local  time  stepping,  within  about  300  iterations,  and  an  argument  could  be  made  that 
convergence  was  reached  in  only  250  iterations,  as  illustrated  by  Figure  4-1. 


Figure  4-1,  Inviscid  Waverider  Convergence  History 


In  comparison,  the  explicit  version  of  the  code  required  over  950  iterations,  with  local 
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time  stepping,  to  reach  convergence  (Figure  4-2).  The  fewer  number  of  iterations 
required  by  the  implicit  method  shows  its  strength. 


Figure  4-2,  Inviscid  Waverider  Convergence  History 


The  CFL  number  shown  in  the  above  figures  is  a  constant  used  in  determining 
the  time  step,  as  shown  by  Equations  (4-2)  through 

=  CFLminj^^A/j  (4-2) 

where 
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The  variable  5  is  the  distance  across  the  cell,  with  the  subscript  indicating  direction.  The 
variable  c  is  the  local  speed  of  sound. 

The  drawback  to  the  implicit  version  is  that  each  iteration  of  the  implicit  method 
takes  considerably  more  computer  time  than  each  iteration  of  the  explicit  method.  The 
implicit  version  of  the  code  took  0.0001  CPU  seconds  per  iteration  j)er  node.  The 
explicit  version  required  approximately  0.000054  CPU  seconds  per  iteration  per  node, 
approximately  twice  as  fast  as  the  implicit  version.  The  explicit  version,  however, 
requires  only  85%  of  the  memory  that  the  implicit  version  does. 

The  shock  capturing  ability  of  both  methods  is  essentially  equivalent.  Contour 
plots  of  density  on  a  Y-Z  plane  near  the  base  plane  (Figure  4-3  and  Figure  4-4)  shows 
a  very  small  difference  between  the  results  of  the  two  methods.  Figure  4-3  and 
Figure  4-4  show  the  results  for  the  Y-Z  plane  just  upstream  of  the  baseplane.  Note  that 
the  shock  is,  effectively,  entirely  c£q}tured  by  the  compression  surface  of  the  waverider. 
There  is  a  slight  expansion,  or  spillage,  at  the  leading  edge  due  to  numerical  truncation 
of  the  waverider  model,  and  the  numerical  error  of  the  solution.  The  two  figures  are 
very  much  alike,  which  is  comforting,  showing  that  the  two  methods  reached  the  same 
solution.  The  implicit  method  seemed  to  do  a  better  job  of  capturing  the  shock,  though 
that,  as  the  different  placement  of  the  contour  line  closest  to  the  waverider  surface,  may 
be  due  to  the  different  convergence  the  two  solutions  went  to.  The  results  shown  in 
Figure  4-4  are  essentially  identical  to  those  shown  in  Figure  4-3.  There  is  a  slight 


difference  in  the  placement  of  three  of  the  contour  lines,  but  that  can  be  attributed  to  the 
different  residuals  of  the  two  cases.  Basically,  the  implicit  and  explicit  methods 
produced  the  same  results.  L/D  differed  in  the  5*  decimal  place,  the  shocks  formed  in 
the  same  place,  and  the  distributions  of  the  flow  variables  were  practically  identical. 
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Figure  4-3,  Implicit  Method,  Inviscid  Waverider 
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4.1.2  Modifled  Waverider 


An  inviscid  computational  solution  for  the  waverider  modified  for  viscous  flow 
gave  results  very  similar  to  the  unmodified  waverider.  A  sparser  grid  than  that  used  for 
the  Navier-Stokes  case  was  used  in  an  effort  to  conserve  computer  resources.  The 
implicit  method  was  used  exclusively  on  the  modified  waverider,  again,  in  an  effort  to 
conserve  computer  resources.  The  L/D  ratio  was  slightly  lower  than  for  the  unmodified 
waverider,  due  mainly  to  a  rise  in  drag.  The  L/D  ratio  was  calculated  to  be  8. 158,  only 
0.4%  less  than  that  calculated  by  Stecklein  (28)  for  the  inviscid  optimized  waverider. 
Both  lift  and  drag  changed,  drag  most  noticeably.  The  change  to  lift  was  practically 
unnoticeable.  A  contour  plot  of  density  at  the  base  plane.  Figure  4-5,  shows  the  shock 
forming  similarly  to  the  unmodified  waverider  model.  Figure  4-3.  The  same  flow 
spillage  seen  in  Figure  4-3  and  Figure  4-4  was  noted  at  the  leading  edges,  although  now 
it  is  caused  more  by  the  perturbation  of  the  flow  by  the  blunt  leading  edge  than  by 
numerical  errors  in  the  waverider  model  and/or  the  solution. 
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Figure  4-5,  Modified  Waverider  Model,  Inviscid  Flow 


The  blunted  nose  produced  a  detached  bow  shock.  This  detached  bow  shock 
produced  very  high  temperatures  in  the  nose  region,  since  the  shock  was,  in  that  region, 
a  normal  shock.  The  shock  quickly  dissipated  on  the  freestream  surface  as  the  flow 
expanded  around  the  nose,  creating  an  expansion  wave  which  intersected  and  canceled 
the  shock,  although  the  shock  affected  the  temperature  and  Mach  number  distribution 
on  the  waverider’s  surface.  On  the  compression  surface,  the  shock  quickly  became  an 
oblique  shock  as  the  flow  expanded  around  the  nose,  but  as  on  the  freestream  surface, 
flow  temperature  and  Mach  number  were  affected,  particularly  in  the  nose  region. 

Figure  4-6  shows  the  temperature  and  Mach  number  variation  over  the  entire 
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freestream  surface  of  the  waverider.  Figure  4-7  shows  a  zoom  of  the  nose  region.  Note 
the  r^id  decrease  in  temperature  away  from  the  nose,  and  the  somewhat  slower,  though 
still  rapid  rise  in  Mach  number  to  the  freestream  value  of  Mach  10.  The  rise  in  Mach 
number  is  due  to  the  combination  of  the  increase  in  flow  velocity  and  the  drop  in 
temperature. 

Figure  4-8  shows  the  Mach  number  and  temperature  distributions  on  the  whole 
of  the  compression  surface,  and  Figure  4-9  shows  a  zoom  of  the  nose  region.  The 
behavior  of  the  flow  on  the  compression  surface  is  similar  to  that  on  the  freestream 
surface,  though  on  the  compression  surface  the  flow  expands  to  match  conditions  behind 
an  oblique  shock.  As  on  the  freestream  surface,  the  Mach  number  increases,  and 
temperature  rapidly  decreases.  On  the  compression  surface,  however,  the  temperature 
is  generally  higher  than  on  the  freestream  surface,  and  the  Mach  number  lower. 
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4.2  Viscous  Flow  Results 

Convergence  for  the  Navier-Stokes  equations  and  the  computational  grid  required 
to  resolve  the  flow  was  quite  slow.  Small  grid  volumes  near  the  leading  edges  and 
nose,  while  providing  good  flow  resolution,  required  very  small  values  for  the  time  step 
in  order  to  keep  from  exceeding  the  allowable  CFL  stability  criteria.  Figure  4-10  shows 
the  convergence  history  for  the  viscous  flow  case. 


Figure  4-10,  Viscous  Flow  Convergence  History 
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The  convergence  went  as  expected,  though  not  as  quickly  as  desired.  The 
residual  at  first  rose  and  then  began  a  short  series  of  cycles  in  which  the  maximum 
residual  decreased.  After  a  significant  decrease  in  residual  another  series  of  cycles 
began  in  which  the  maximum  residual  reached  in  the  cycle  decreased.  The  CFL 
stability  criteria  was  adjusted  automatically  by  the  CFD  code  as  necessary.  At  first  the 
CFL  number  was  quite  low,  on  the  order  of  0.09,  but  then  rose  rapidly.  The  CFL 
number  did  not  remain  steady  at  the  maximum  allowed  until  approximately  the  800* 
iteration.  The  root  mean  square  surface  pressure,  or  pressure  residual,  behaved  as 
expected,  rising  rapidly  and  then  leveling  off.  A  slight  decrease  was  noted,  and,  after 
a  brief  period  of  little  to  no  change,  the  pressure  residual  began  to  rise  slowly.  A 
constant  pressure  residual  is  a  good  indication  that  the  solution  is  converged,  as  the 
pressure  on  the  waverider  surface  is  not  changing  significantly,  and  thus  the  changes  to 
the  flow  structure  with  each  iteration  are  very  minor. 

As  expected,  when  viscous  terms  are  included  in  the  equations,  L/D  drops  due 
to  the  addition  of  viscous  drag.  The  L/D  ratio  for  the  waverider  in  viscous  flow  is  5.74, 
a  30%  drop  from  the  inviscid  case.  This  value  of  L/D  is  similar  to  that  calculated  in 
other  solutions  of  viscous  flow  over  a  waverider.  Another  primary  contribution  to  drag 
is  wave  drag,  which  is  caused  by  the  shock  waves.  In  this  case  a  weak  shock  formed 
over  the  freestream  surface  of  the  waverider,  as  illustrated  by  the  lower  part  of 
Figure  4-11.  This  shock  was  formed  in  part  by  the  bow  shock  produced  by  the  blunt 
leading  edges,  and  in  part  by  the  formation  of  the  boundary  layer.  The  shock  on  the 
compression  surface  was  also  affected  by  the  formation  of  the  boundary  layer.  The 
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boundary  layer  on  the  compression  surface  was  thinner  than  on  the  freestream  surface, 
but  had  the  effect  of  increasing  the  effective  flow  turning  angle,  increasing  shock 
strength  slightly.  The  stronger  shocks,  caused  by  the  interaction  of  the  boundary  layer 
and  shock  waves,  contributed  to  the  drop  in  L/D  as  compared  to  the  inviscid  flow  case 
by  increasing  the  wave  drag.  Thus  viscosity  contributed  to  drag  in  two  ways.  It 
contributed  directly  and  by  increasing  the  strength  of  the  shocks,  and  thus  the  wave 
drag.  The  computed  L/D  ratio  compared  favorably  with  results  for  similar  vehicles. 
The  computed  L/D  was  greater  than  that  determined  by  Muller  (22),  who  arrived  at  a 
value  of  4.18  for  a  waverider  flying  at  Mach  5.5,  though  substantially  less  than  those 
reported  by  Chang  (6)  and  Takashima  (31).  Compared  to  Rasumussen’s  (26)  analytic 
work,  which  gave  results  veiy  like  those  reported  by  Muller,  the  computed  value  was 
high,  though  not  unreasonably  so.  Comparing  the  result  of  5.74  obtained  in  this 
research  to  a  plot  of  L/D  such  as  that  in  Anderson  (5),  this  result  can  be  seen  to  be  quite 
near  the  predicted  value,  and  the  results  obtained  by  others. 

Figure  4-11  illustrates  the  pressure  field  at  the  waverider’s  line  of  symmetry. 
Figure  4-11  shows  the  oblique  shock  that  forms  below  the  compression  surface  (though, 
since  the  waverider  is  depicted  herein  upside  down,  the  shock  appears  to  be  above  the 
waverider),  and  the  shock  and  expansion  that  form  above  the  freestream  surface. 
Figure  4-12  shows  density  contours  for  the  same  region  as  that  depicted  in  Figure  4-11. 
The  shock  on  the  compression  side  is  clearly  depicted,  as  is  the  much  weaker  shock  on 
the  freestream  side.  Figure  4-13  is  a  zoom  of  the  density  contours  in  the  nose  region. 
This  figure  shows  the  shock  on  the  freestream  side  more  clearly  than  does  Figure  4-12. 
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Figure  4-14  shows  the  Mach  number  contours  in  the  nose  region.  The  expansion  can 
be  seen  here  as,  just  above  the  freestream  surface,  the  flow  accelerates.  Note  also  the 
rapid  increase  in  Mach  number  above  the  surface,  indicating  how  thin  the  boundary 
layer  is  in  the  nose  region.  Figure  4-15  is  a  further  zoom  of  the  Mach  contours  in  the 
region  of  the  waverider’s  nose.  This  figure  shows  even  better  than  the  previous  one 
how  thin  the  boundary  layer  is  in  this  region.  It  also  shows  the  shock  structure  in  the 
nose  region  quite  clearly. 
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Figure  4-12,  Centerline  Density  Contours 
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The  velocity  vector  plots.  Figure  4-16  through  Figure  4-18,  show  the  steep 
velocity  gradients  in  the  boundary  layer  quite  well.  The  boundary  layer  grew  more 
quickly  on  the  freestream  surface  than  on  the  compression  surface  due  to  the  lower 
pressures.  Comparing  Figure  4-17  with  Figure  4-18  it  can  be  seen  that,  at  the 
waverider’s  trailing  edge,  the  boundary  layer  on  the  freestream  surface  is  almost  twice 
as  thick  as  on  the  compression  surface. 

The  temperature,  as  shown  in  the  plots  of  temperature  and  density  vs.  Y 
(Figure  4-19  through  Figure  4-22),  changes  as  predicted  by  theory  for  a  cold  wall; 
climbing  deeper  into  the  boundary  layer  to  a  maximum  just  above  the  surface,  and  then 
falling  off  sharply  to  the  surface  temperature.  In  the  plots  of  temperature  presented  here, 
the  temperature  did  not  reach  the  surface  temperature  as  the  data  was  calculated  at  the 
first  cell  center,  which  was  positioned  just  above  the  surface  of  the  waverider. 
Interpolation  between  the  value  at  the  first  point  above  the  waverider  surface  and  the 
ghost  point  just  below  the  surface  returns  the  value  for  temperature  imposed  as  a 
boundary  condition.  Figure  4-19  shows  more  the  effects  of  the  shock  and  expansion 
than  the  boundary  layer,  due  to  the  boundary  layer’s  thinness  near  the  nose.  Figure  4-20 
also  shows  more  the  effects  of  the  expansion  than  the  boundary  layer.  The  sudden  turn 
in  the  plot  at  Y  =  0.025  feet  is  due  to  the  thinness  of  the  boundary  layer.  The  boundary 
layer  is  entirely  missed,  and  the  flow  almost  jumps  to  the  imposed  boundary  conditions. 

The  plots  of  flow  variables  in  the  boundary  layer  near  the  nose  region  show  how 
thin  the  boundary  layer  is  there.  Also,  because  the  shock  angle  is  so  small,  the  plots 
extend  almost  to  the  shock.  A  much  finer  grid  would  be  required  to  resolve  the 
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boundary  layer  at  all  near  the  waverider  nose. 

The  plots  of  temperature  and  density  near  the  trailing  edge.  Figure  4-21  and 
Figure  4-22,  show  the  expected  behavior,  described  above,  very  well.  Both  plots  show 
temperature  rising  sharply  just  above  the  waverider  surface  to  a  maximum  that  is 
approximately  at  the  center  of  the  boundary  layer,  and  then  falling  off  to  the  temperature 
at  the  edge  of  the  boundary  layer.  Density  mirrored  the  behavior  of  temperature, 
indicating  an  almost  constant  pressure  through  the  boundary  layer. 

The  plots  of  velocity  and  Mach  number.  Figure  4-23  through  Figure  4-26,  mirror 
the  behavior  seen  in  the  vector  plots.  Figure  4-17  and  Figure  4-18.  As  the  wall  is 
approached  the  velocity  drops  rapidly  to  zero.  As  explained  above,  the  values  shown 
don’t  actually  reach  zero  velocity,  as  the  values  are  not  calculated  on  the  surface  of  the 
waverider,  but  rather,  just  above  it.  Interpolation  between  the  ghost  point  and  the  first 
point  above  the  waverider  surface  shows  that  both  the  velocity  and  Mach  number  are 
zero  at  the  surface.  The  plots  of  velocity  and  Mach  number  near  the  nose.  Figure  4-24 
and  Figure  4-23  show  the  effects  more  of  the  flow  expansion  in  that  region  than  the 
boundary  layer.  The  plot  of  the  velocity  profile  above  the  compression  surface  shown 
in  Figure  4-23  looks  like  a  boundary  layer  profile,  but  is  actually  caused  by  the  flow 
expansion  in  the  region.  The  plot  of  the  velocity  profile  above  the  freestream  surface. 
Figure  4-24,  shows  a  very  thin  boundary  layer.  Like  Figure  4-20,  the  values  "jump", 
that  is  the  slope  is  discontinuous,  at  the  points  just  above  the  waverider  surface.  This 
is  due  to  inadequate  resolution  of  the  flow  in  the  boundary  layer,  and  the  boundary 
conditions  being  imposed,  just  as  explained  above  for  the  plots  of  temperature  and 
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Figure  4-16,  Centerline  Velocity  Vector  Plot,  Nose  Region 
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Figure  4-17,  Centerline  Velocity  Vector  Plot,  Trailing  Edge,  Compression  Surface 
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Figure  4-18,  Centerline  Velocity  Vector  Plot,  Trailing  Edge,  Freestream  Surface 
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Figure  4-19,  Compression  Surface  Centerline  Boundary  Layer  Data,  Nose  Region 
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Figure  4-21,  Compression  Surface  Centerline  Boundary  Layer  Data,  Trailing  Edge 
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Mach  Number  &  Velocity  vs  Y 
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Figure  4-23,  Compression  Surface  Centerline  Boundary  Layer  Data,  Nose  Region 
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Figure  4-25,  Compression  Surface  Centerline  Boundary  Layer  Data,  Trailing  Edge 
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Figure  4-26,  Freestream  Surface  Centerline  Boundary  Layer  Data,  Trailing  Edge 
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Heat  transfer  from  the  flow  to  the  waverider  was  greatest  at  the  nose  and  leading 
edges,  as  was  the  skin  friction  coefficient.  Figure  4-27  and  Figure  4-28  illustrate  the 
heat  transfer  and  skin  friction  compression  surface  and  Figure  4-29  and  Figure  4-30 
show  these  values  on  the  freestream  surface.  The  variation  of  heat  transfer  and  skin 
friction  was  greatest  in  the  spanwise  direction,  as  can  be  seen  in  the  following  figures. 
The  important  thing  to  note  in  these  figures  is,  as  mentioned  above,  that  the  skin  friction 
and  heat  transfer  are  greatest  at  the  nose  and  leading  edges. 
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Figure  4-27,  Compression  Surface  Stanton  Number  Distribution 
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Figure  4-28,  Compression  Surface  Skin  Friction  Coefficient  Distribution 
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Figure  4-3 1  illustrates  the  variation  of  pressure,  temperature,  and  velocity  on  the 
stagnation  stream  line.  The  shock  standoff  from  the  nose  is  approximately  0.02  inches, 
a  very  small  distance,  indicating  that  the  stagnation  region  will  probably  comprise  a 
significant  portion  of  this  distance.  As  shown  by  Figure  4-31,  the  flow  was  not  well 
resolved  in  the  stagnation  region.  A  much  finer  grid  would  be  required  to  accurately 
resolve  the  flow  in  the  stagnation  region. 

The  figure  does  show,  by  inference,  the  extreme  flow  gradients  present  in  the 
stagnation  region.  The  flow  will  have  to  slow  from  a  little  more  than  5000  feet  per 
second  to  zero  velocity  within  .025  feet.  This  will  cause  a  corresponding  rise  in 
temperature  as  the  kinetic  energy  is  transformed  to  internal  energy.  Temperature  will 
show  extremely  steep  gradients,  as  it  not  only  reaches  a  maximum,  but  then  has  to  drop 
to  match  the  wall  temperature.  These  extreme  gradients,  combined  with  too  coarse  of 
a  grid,  make  the  numerical  solution  in  the  stagnation  region  to  be  of  questionable  value 
and  accuracy. 
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A  Y-Z  cross  section  of  the  flow  near  the  baseplane.  Figure  4-32,  looked  quite 
similar  to  those  of  the  inviscid  case,  as  shown  in  Figure  4-3  and  Figure  4-4.  The  main 
differences  between  the  viscous  and  inviscid  cases  included  the  boundary  layer,  a 
compression  above  the  freestream  surface  due  to  the  boundary  layer,  and  the  shock 
detaching  from  the  leading  edges  of  the  waverider  to  stand  off  a  short  distance.  The 
baseplane  showed  an  interesting  phenomena;  a  lambda  shaped  shock,  shown  in 
Figure  4-33.  This  phenomenon  was  caused  by  a  flow  expansion  around  the  waverider’ s 
leading  edge  meeting  the  high  pressure  region  on  the  compression  surface.  The  strong 
compression  from  the  low  pressure  just  inboard  of  the  leading  edge  to  the  pressure  on 
the  main  part  of  the  compression  surface  caused  the  shock  to  form.  This  shock  structure 
did  not  exist  two  planes  upstream  of  the  baseplane;  a  single  detached  shock  appears,  as 
pictured  in  Figure  4-34,  which  is  two  planes  upstream  of  the  baseplane. 
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Figure  4-3  ♦,  Density  Contours,  Viscous  Flow,  50“*  plane  from  Nose,  Zoom  of  Leading 
Edge 


93 


V.  THEORETICAL  AND  NUMERICAL  ANALYSIS 


5.1  Boundary  Layer  Equations 

The  compressible  boundary  layer  equations  are  derived  directly  from  the  Navier- 
Stokes  equations  through  order  of  magnitude  analysis.  Rasmussen  (27)  derives  the  axi- 
symmetric,  steady  state,  compressible  boundary  layer  equations  in  some  detail,  with  the 
final  result  b  .ing; 

Continuity:  =  0  (5-1) 

dx  dy 

Where  m  is  0  for  planar  two-dimensional  flow  and  1  for  axisymmetric  flow  The  only 
difference  between  the  planar  equations  and  the  axisymmetric  equations  is  the  r"  term 
in  the  continuity  equation. 

X-Momentum:  +  pv—  =  -  -El  +  —  [■  (5-2) 

dx  dy  dx  dy J 

Y -Momentum:  EE.  =  0  (5-3) 

dy 
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and  the  integral  enthalpy  equation  is 


(5-5) 
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(5-6) 


Where  and  are  total  enthalpy  at  the  wall  and  boundary  layer  edge,  respectively,  and 
St  is  Stanton  number.  The  terms  8*,  displacement  thickness,  and  6,  momentum 
thickness,  are  given  by 


and  St  is  defined  as 
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Equations  (5-5)  and  (5-6)  can  be  related  through  Reynold’s  analogy.  If  Prandtl 
number,  Pr,  is  equal  to  one,  Reynolds  analogy  states  that  the  left  hand  sides  of  equations 
(5-5)  and  (5-6)  are  equal,  so  that  heat  transfer  is  directly  related  to  skin  friction.  If  Pr 
is  not  equal  to  one,  Reynold’s  analogy  can  be  written  as  (3) 
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(5-9) 
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2Pr^ 

For  analysis  of  the  stagnation  region  near  the  nose,  a  similarity  solution  is  useful.  Such 
an  analysis  can  be  accomplished  through  the  use  of  the  Levy-Lees  transformation. 
Anderson  (3)  provides  a  transformation  to  the  Levy-Lees  coordinates,  ^  and  r\: 

X 

^  =  JP.". 
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T1  = 


(5-11) 


Anderson  (3)  derives  the  similarity  form  of  the  boundary  layer  momentum  and 
energy  equations: 

X-Momentum: 


4-icJ)  +//  = 
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Y-Momentum: 
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(5-12) 


and  Energy: 


^  =0 
dr] 


(5-13) 
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where 


*  /g  =  2£  -  gK  *  -  A 

^Fi-*J  -'*  u,cK,  *a^  pA,  0^  A/ 


g(4.n)  =  -r- 
h 


(5-15) 


^(y(5,n))  =  /  =  - 


(5-16) 


(5-17) 


A  computer  code  (18)  which  numerically  implements  the  similarity  solution, 
described  in  equations  (5-12)  through  (5-14),  for  the  stagnation  region  of  hypersonic 
flow  around  a  blunt  ended  cylinder  was  used  to  compare  the  flow  in  the  stagnation 
region  to  an  analytic  solution.  The  following  two  figures  are  plots  of  data  produced  by 
the  similarity  program.  Figure  5-1  graphically  displays  the  extreme  gradients  typical  of 
hypersonic  flow.  The  flow  is  slowed  from  1735  feet/second  to  0  feet/second  in  the 
space  of  less  than  half  an  inch.  The  edge  velocity  is  assumed,  here,  to  be  equal  to  that 
just  behind  the  normal  shock. 

The  similarity  solution  does  not  match  particularly  well  with  the  data  taken  from 
the  computational  solution  (Figure  4-31).  This  is  primarily  due  to  poor  flow  resolution 
in  the  stagnation  region.  A  finer  computational  grid  in  the  stagnation  region  should 


Figure  5-1,  Stagnation  Line  Plot  of  Velocity  and  g 

correct  the  difficulty.  Hayes  and  Probstein  (19)  predict  that  the  ratio  of  shock  stand-off 
distance  to  nose  radius  as  0.104,  for  M«,=10  and  y  =  1.4.  This  matches  well  with  vdiat 
observed  in  the  computational  solution,  but  also  means  that  due  to  the  very  small 
shock  stand-off  distance,  the  stagnation  region  will  occupy  a  large  portion  of  the  shock 
l^er. 

5.1.1  Hypersonic  Viscous  Interaction 

Anderson  (3)  states  that  the  boundary  layer  thickness  in  a  viscous,  hypersonic 
flow  over  a  flat  plate  is  proportional  to  Mach  number  squared  divided  by  the  square  root 
of  Reynolds  number.  He  gives  the  relation: 
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(5-18) 


S  oc 


The  waverider,  locally,  resembles  a  flat  plate,  particularly  on  the  freestream 
surface  which  is  parallel  to  the  freestream.  The  main  difficulty  in  applying  equation 
(5-18)  to  the  compression  surface  is  its  the  angle  of  attack.  However,  since  equation 
(5-18)  is  not  an  equality  but  a  proportionality,  it  should  still  hold,  though  the  constant 
of  proportionality  will  be  different  than  for  a  flat  plate  at  a  zero  degree  angle  of  attack. 

The  rapid  boundary  layer  growth  and  the  resulting  thickness  of  the  boundary 
layer  are  the  basic  drivers  of  hypersonic  viscous  interaction.  Anderson  (3)  goes  on  to 
derive  the  equation 


.^1.  =  1  *  <I,X  (S-l») 

p. 

for  strong  viscous  interactions,  where  a,  is  a  constant.  Note  that  equation  (5-19) 
indicates  that  the  pressure  ratio,  p,/pm  is  dependent  only  on  x,  and  that  the  pressure 
ratio  is  linear  with  x- 

For  a  weak  hypersonic  viscous  interaction  Anderson  (3)  derives 

IL‘  \  *  by,  (5-20) 

where  b,  is  a  constant.  Again,  the  pressure  ratio  varies  linearly  with  x-  Anderson  (3) 
provides  numerical  values  for  the  constants  ay  and  b,  for  hypersonic  flow  over  both  an 
adiabatic  wall  and  a  cold  wall  flat  plate.  Equations  (5-19)  and  (5-20)  become: 
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-^  =  1  +  0.5x 


{strong  interaction) 


(5-21) 


_  =  1  +  0.078x  {weak  interaction) 

Po. 

for  the  cold  wall  case.  With  the  viscous  interaction  parameter  x  >  3.0,  strong  interaction 
is  generally  taken  to  occur,  and  weak  interaction  for  x  <  3.0  (3).  Note  that  the  strong 
interaction  case  varies  much  more  rapidly  with  %  than  does  the  weak  interaction  case 
because  the  constant  multiplying  x  is  an  order  of  m^nitude  larger  for  the  strong 
interaction  case  than  for  the  weak  interaction  case. 

The  definition  of  %  in  the  above  equations  is. 

(S-22) 


where  C  is  given  by 


C  =  (5-23) 

P.l*. 


5.1.2  Application  of  the  Boundary  Layer  Equations 

The  above  equations  were  applied  to  the  waverider  flying  at  Mach  10  and  a  = 
0  at  an  altitude  of  100,000  feet.  Figure  5-2  shows  the  estimated  variation  of  x  with 
distance  from  the  waverider’s  nose.  For  purposes  of  calculating  x  in  Figure  5-2  the 
following  assumptions:  u,  =  =  8910  feet/second,  T,  =  1469  ®Rankine,  and  C  =  1.25 

were  made.  Sutherland’s  formula  was  used  to  calculate  and  Re  was  calculated  using 
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the  value  of  the  variables  at  the  boundary  layer’s  edge. 


Figure  5-2,  Chi  vs  X 

Applying  %  as  shown  in  Figure  5-2  to  equation  (5-21),  the  local  pressure  ratio 
is  calculated.  The  pressure  ratio  distribution  over  x,  calculated  both  from  the 
computational  data  and  from  the  viscous  interaction  equations  given  above,  is  shown  in 
Figure  5-3  (compression  surface)  and  in  Figure  5-4  (freestream  surface). 
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Figure  5-3,  Compression  Surface  Centerline  Pressure  Ratio 
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Figure  5-4,  Freestream  Centerline  Pressure  Ratio, 
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Appropriate  flow  variables,  particularly  were  extracted  from  the  computational 
solution  for  use  in  the  calculation  of  x- 

As  can  be  seen  in  both  flgures,  the  match  is  not  particularly  good.  The  initial 
trend  of  a  steep  decrease  in  pressure  ratio  is  correct,  though  the  values  don’t  match.  The 
numerical  results  are  initially  greater  than  the  theoretic  prediction.  The  steep  pressure 
rise  behind  the  normal  shock  accounts  for  the  displacement  of  the  computational 
pressure  ratio  above  the  predicted  ratio.  The  sudden  turn  and  then  rise  in  the 
computationally  calculated  pressure  ratio  is  due  to  the  flow  expanding  past  freestream 
conditions  around  the  blunt  nose,  and  then  trying  to  recover  to  the  freestream  pressure, 
as  described  earlier,  in  Chapter  4. 


5.2  Skin  Friction  and  Heat  Transfer 

Incompressible  flow  over  a  flat  plate  produces  the  following  skin  friction  and 
heat  transfer  coefficients  (36)  derived  from  similarity  solutions; 


0.664 


(5-24) 


and 


St 


0.332 


(5-25) 


The  parameters  Re  and  Pr  are  based  on  the  flow  properties  at  the  edge  of  the  boundary 
layer. 
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The  reference  temperature  method  uses  the  exact  same  form  for  the  skin  friction 
and  heat  transfer  coefficients  as  is  found  in  equations  (5-24)  and  (5-25),  except  that  Re 
and  Pr  are  now  calculated  using  a  reference  temperature,  T*,  given  by  Anderson  (3)  as: 


+  0.032M/  +  0.58 


V  «  J 


(5-26) 


Thus,  Re  and  Pr  can  be  written  as: 


(5-27) 


The  values  in  equation  (5-27)  are  then  substituted  directly  into  equations  (5-24) 
and  (5-25).  The  reference  values  used  in  equation  (5-27)  were  taken  from  the  just 
behind  the  shock  at  about  the  midpoint  of  the  waverider,  to  ensure  that  the  reference 
values  used  were  indeed  representative  of  the  conditions  at  the  edge  of  the  boundary 
layer. 


The  plot  of  skin  friction  coefficient  versus  x  (feet)  in  Figure  5-5  shows  results 
very  similar  to  those  returned  by  the  computational  solution,  shown  in  Figure  5-6,  The 
agreement  between  the  two  solutions  is  good.  Both  figures  show  precipitous  drop  in 
both  heat  transfer  and  skin  fHction  as  distance  from  the  nose  increases.  Both  plots 
asymptotically  approach  a  value  of  0.00004  for  skin  friction  coefficient  and  .00008  for 
Stanton  number.  Agreement  with  Reynold’s  analogy  in  the  numerical  solution  is  also 
good,  Stanton  number  being  a  little  less  than  twice  the  value  of  the  skin  faction 
coefficient. 
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Figure  5-5,  Cf  versus  x 


Figure  5-6,  St  &  G)inpression  Surface  Centerline 
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5.3  Lift  and  Drag 

The  lift  and  drag  forces  are  calculated  by  integration  of  the  pressure  and  shear 
stresses  over  the  body.  White  (36)  gives  the  control  volume  equation; 


fff^p 

JjdA  JJ^dm 


dVol 


(5-28) 


where  F,  is  the  surface  force  and  F*  is  the  body  force.  The  total  force,  F,  is  then 
broken  into  lift  and  drag  components  by  simply  taking  lift  as  the  component  of  F  in  the 
negative  Y  direction  (see  Figure  1-2)  and  drag  as  the  component  in  the  direction  of  the 
freestream  flow.  A  program,  written  by  Gaitonde,  which  analyzes  the  data  produced  by 
the  CFD  solver  further  breaks  the  lift  and  dr£^  components  into  viscous  and  pressure 
components.  Table  S-1  lists  the  forces  for  the  viscous  and  inviscid  cases  normalized  to 
the  force  in  the  X  direction. 


Table  5-1,  Force  Data 


Parameter 

Viscous  Flow  Value 

Inviscid  Flow  Value 

Total  Lift,  lbs 

223198 

195805 

Total  Drag,  lbs 

38853 

23999 

Pressure  Lift,  lbs 

223639 

195805 

Pressure  Drag,  lbs 

27708 

23999 

Viscous  Lift,  lbs 

-441 

0 

Viscous  Drag,  lbs 

11145 

0 

149.938 

142.593 

P«,  psf 

23.276 

23.276 

Re, 

193,043,000 

L/D 

5.745 

8.159 
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The  flight  conditions  outlined  in  Table  4-1  were  used  in  the  calculation  of  both  the 
viscous  and  inviscid  cases. 

As  can  be  seen,  the  viscous  terms  contribute  significantly  to  the  drag  on  the 
vehicle,  almost  30%  of  the  drag  being  due  to  viscosity.  The  viscous  effects  increase 
lift  by  12%  as  compared  to  the  inviscid  case.  The  increase  in  lift  is  due  to  the  slightly 
greater  pressure  on  the  compression  surface,  caused  by  the  hypersonic-boundary  layer 
interaction  present  in  the  viscous  flow  solution. 
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VI.  CONCLUSIONS  AND  RECOMMENDATIONS 


6.1  Conclusions 

Mach  10  Navier-Stokes  flow  over  a  conically  derived  hypersonic  waverider  was 
computationally  solved  on  the  OSC  CRAY  Y-MP  using  WL/FIM’s  implicit  flux 
splitting  code.  The  code  performed  well,  its  finite  volume  formulation  enhancing  the 
practical  stability,  despite  some  very  skewed  cell  geometries. 

Returning  to  the  objectives  stated  in  chapter  I,  the  following  conclusions  were 
reached: 

Blunting  the  nose  and  leading  edges  had  little  effect  on  results  such  as  lift  to 
drag  ratio,  though  it  did  significantly  effect  the  flow  structure,  particularly  near  the  nose 
of  the  vehicle. 

Generating  the  three  dimensional  grid  was  a  major  effort,  though  made 
considerably  simpler  by  the  GRIDGEN  package.  The  use  of  the  0-grid  allowed  for 
good  resolution  near  the  waverider’s  leading  edges,  allowing  the  capturing  of  flow 
expansion  and  shock  standoff. 

The  three  dimensional  implicit  flux  splitting  algorithm  performed  extremely  well. 
The  finite  volume  ^proach  the  code  used  made  for  an  extremely  robust  algorithm,  able 
to  handle  regions  with  highly  skewed  cells,  and  cells  with  very  small  volumes.  The 
code  produced  results  predicted  by  hypersonic  theory,  with  deviations  attributable  to 
numerical  error  in  either  the  waverider  model  or  in  the  numerical  solution  procedure. 
The  implicit  version  of  the  code  required  fewer  iterations  to  reach  convergence  than  did 
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the  explicit  version.  The  explicit  version,  however,  required  less  computation  time  per 
iteration  and  less  memory  overhead.  The  tradeoffs  between  the  versions  of  the  code 
favored  the  implicit  version.  While  each  iteration  for  the  implicit  version  required 
^proximately  twice  as  much  computer  time,  the  implicit  version  required  approximately 
one  quarter  the  number  of  iterations.  If  memory  becomes  a  limiting  factor,  however, 
the  explicit  version  of  the  code  has  a  definite  advantage. 

The  blunted  nose  and  leading  edges  led  to  a  small,  but  strong,  normal  shock 
slightly  displaced  from  the  waverider.  The  pressure  rise  through  the  normal  shock,  and 
the  expansion  around  the  shoulder  of  the  waverider  caused  the  flow  to  differ 
significantly  from  the  results  predicted  by  the  viscous  interaction  equations  given  by 
Anderson  (3),  making  the  comparison  of  little  value.  The  structure  of  the  shock  waves 
was  simple  and  held  only  one  minor  surprise,  the  X  shaped  structure  at  the  leading  edge 
on  the  baseplane,  formed  by  an  expansion  around  the  leading  edge  meeting  the  high 
pressure  on  the  compression  surface.  Elsewhere,  the  shock  was  detached  from  the 
waverider  by  a  small  distance. 

Comparison  of  the  stagnation  region  near  the  nose  to  a  similarity  solution  for 
flow  over  a  similar  geometry  yielded  good  results.  The  trends  matched  well,  and  the 
values  computed  were  close  to  those  predicted.  Differences  can  easily  be  attributed  to 
differences  in  the  geometry,  and  highlight  the  need  for  a  finer  grid  in  the  stagnation 
region. 

The  computed  L/D  ratio  compared  favorably  with  results  for  similar  vehicles. 
Comparing  the  result  of  5.75  obtained  in  this  research  to  a  plot  of  L/D  such  as  that  in 
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Anderson  (5),  this  result  can  be  seen  to  be  quite  near  the  predicted  value,  and  the  results 
obtained  by  others. 

6.2  Recommendations 

This  thesis,  along  with  Stecklein’s  (28),  will  provide  an  excellent  basis  from 
which  to  continue  waverider  research.  The  following  recommendations  are  the  result 
of  careful  consideration  of  the  research  and  conclusions  reached  in  the  course  of  this 
thesis. 

The  use  of  a  taut  cubic  spline,  or  radius  method  to  blunt  the  nose  and  leading  edge  is 
recommended.  Such  a  method  will  provide  a  smoother  model  geometry  at  the  nose  and 
leading  edges  than  that  used  here.  ReHning  the  computational  grid  is  necessary, 
particularly  in  the  stagnation  region  near  the  waverider’s  nose.  Some  particular  goals 
to  pursue  include  reducing  the  number  of  cells  in  the  grid  without  losing  the  flow 
resolution,  and  reducing  the  cell  skewness,  particularly  near  the  leading  edges. 

This  research  neglected  the  real  gas  effects,  such  as  dissociation,  chemical  reaction,  and 
ionization.  Murthy  (23),  and  or  Anderson  (3),  describe  these  effects  the  influence  they 
have  on  the  flow,  but  the  most  noticeable  effects  are  a  reduction  in  the  boundary  layer 
thickness  and  in  aerodynamic  heating.  Turbulence  should  be  accounted  for  in  further 
studies.  The  Reynolds  number  indicates  that  turbulence  could  play  a  major  role  in  the 
flow  over  the  waverider. 

Flight  at  off-design  conditions  should  be  examined.  Waveriders  are  very 
sensitive  to  off-design  conditions,  and  pay  heavy  penalties  in  terms  of  efficiency  when 
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flying  at  other  than  their  design  point.  Particularly,  waverider  performance  over  a  range 
of  Mach  numbers  from  0  to  somewhat  above  the  design  Mach  number  should  be 
investigated,  as  should  performance  for  a  range  of  angles  of  attack  and  yaw.  Waverider 
performance  over  a  range  of  altitudes,  peih^s  a  reentry  trajectory,  should  be 
investigated. 
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APPENDIX  A:  DERIVATION  OF  THE  WAVEREDER  EQUATIONS 


The  waverider  equations  are  well  known,  as  is  their  derivation.  However,  most 
derivations  leave  large  steps  to  the  reader’s  imagination  and  frustration.  Therefor,  the 
equations  are  derived  here  in  some  detail.  The  derivations  came  primarily  from 
Rasmussen  (25  and  26)  and  Stecklein  (28). 

The  conically  derived  waverider  equations  come  from  the  Taylor-MacColl 
equation  and  it’s  hypersonic  small  disturbance  theory  form  and  results. 
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Equation  (A-1)  is  the  Taylor-MacColl  equation,  which  produces  the  HSDT  form; 


dV 

- l+cote — i+2V  =0 

ae^  ae 


(A-2) 


Equation  (A-2)  has  analytic  solutions  for  u  and  v  (Vp  and  Vg  respectively), 
which,  using  2nd  order  accurate  T^lor  series  expansion,  are; 
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The  streamlines  between  the  shock  and  the  cone  (i.e.  within  the  shock  layer)  are 
given  by 


V^ds  =  0 


(A-4) 


where  ds  is  differential  length  along  a  streamline,  and  V  is  the  velocity  vector.  The 
equation  says  that  the  velocity  vector  is  parallel  to  the  streamline. 

Expanding  (A-4)  results  in  only  one  component,  that  in  the  direction; 


det 


^ 

«  V  0 
dr  rdQ  0 


={urdQ  -vdr)i^  =  0 


(A-5) 


Rearranging  (A-S)  by  getting  rid  of  the  unit  vector  and  separating  the  variables 


you  get: 


^=ldQ  (A-6) 

r  V 

Further  assuming  that  u«V«,  and  is  furthermore  constant  within  the  shock  layer, 
(A-6)  can  be  integrated  (after  substituting  for  v  from  (A-3))  from  the  shock  to  a  point 
of  interest: 
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(A-7) 


'cdr  ^  _  V  QdQ 

l~r  ^02-6^ 

(A-7)  can  be  solved  to  become: 

Inr  -  Inr^  =Iln|p"  -Ilnje"  -8^1  (A-8) 

which  can  be  rearranged  by  raising  both  sides  to  the  e  to; 

(A-9) 

Now,  the  distance  from  the  centerline  of  the  cone  to  any  given  point  is  rsi>i0  or, 
for  small  angles,  rQ.  The  (arbitrary)  cylindrical  freestream  surface  can  be  expressed  by: 

rsin0=/<j>)  (A-10) 

where ^<t>)  is  arbitrary.  Applying  the  small  angle  approximation,  the  freestream  surface 
which  intersects  the  shock  at  r,  is  given  by: 

/0=r,((|>)P  (A-11) 

In  the  base  plane  we  non-dimensionalize  by  lb,  the  radius  of  the  generating  cone 
in  the  base  plane.  This  gives  =  Q^b.  If  (A-11)  is  evaluated  in  the  base  plane  it 
yields: 

Defining  a  as  p/8  (A-12)  can  be  rearranged  to  get; 
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(A-13) 


CT 

Now,  R^((|>)=10d,((j))/15=0d,((t>)/5,  or  e^=5Rd,((l>).  Evaluating  (A-9)  in  the 
baseplane,  with  (A-13)  substituted  in  for  r,(<|>),  yields: 


eL-6^ 

p^-6^ 


(A-14) 


Squaring  both  sides,  rearranging  slightly,  and  substituting  in  for  0^*: 


(A-IS) 


A  little  algebra  gives  the  final  result: 


= 1  ^ 


a^-l 


(A-16) 


Note  that  is  an  arbitrary  function,  usually  (and  in  this  thesis)  taken  from  the 
transformation  of  the  polynomial: 


X=R^+AY^+BY*+CY^ 


(A-17) 
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APPENDIX  B:  DIVIDED  DIFFERENCES 


[]pi 

[,]pi 

[ , ,  ]Pi 

[ . ,  ,]Pi 

■ 

g(xi) 

Si 

Hi 

g(^i) 

(h.^i+i]g-S|)/AXi 

HZ 

(Si+,+Si-s[Xi,Xi^,]g)/(AXi)^ 

g(Vi) 

(sKt-[ti,tH,]g)/AXi 

Si.i 

Vi 

g(Vi) 

■■ 

■■■■ 

Table  B-1,  Divided  DifTerence  Table 


The  square  brackets,  [  ],  indicate  a  divided  difference  of  order  i-1.  The  first 
divided  difference  is  given  in  Eqn  (B-2)  below. 


fx  -/o 


X  ^  X 


The  second  divided  difference  incorporates  the  first: 


(B-2) 


(B-2) 


Similarly  for  the  third  divided  difference. 

The  values  given  in  Table  B-1  are  used  to  derive  the  coefficients  in  Eqn  (2-1 1). 
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APPENDIX  C:  PROGRAM  WAVERIDER 
c 

c  5  Sep  92: 

c  This  attempt  will  now  proceed  as  follows: 
c  Both  freestream  and  compression  surfaces  will  be  extended  as  y=y(z). 
c  The  slopes  of  each  will  be  clamped  at  -3  for  the  compression  and  +3 
c  for  the  freestream  surfaces,  at  the  boundary  point  x5,y5.  However, 
c  the  boundary  point  value  will  not  actually  be  written  to  the  file, 
c  thus  blunting  the  sharp  comer  that  would  otherwise  result, 
c  The  other  major  change  is  that  the  first  plane,  the  i=l  plane, 
c  will  not  be  printed  to  the  output  file.  Instead,  I3G- VIRGO  will  be 
c  used  to  generate  the  cubic  spline  on  the  nose.  This  will  hopefully 
c  eliminate  singularity  problems  near  the  nose, 
c  Note  also  that  the  slopes  may  be  input,  as  can  the  separation 
c  distance  and  the  number  of  points  to  use  in  the  spline.  Good  results 
c  were  obtained  using  3cm  in  y,  1cm  in  z,  slopes  of  +-4,  and  15  points, 
c 

c  23  Sep  92: 

c  The  program  currently  uses  linear  interpolation  to  build  the 
c  nose  of  tilie  waverider. 
c 

PROGRAM  WVRIDR 

C  PRODUCES  MULTIPLE  2-D  CROSSPLANE  SURFACES  OF  A 
C  PARABOLIC  HYPERSONIC  WAVERIDER 

IMPLICIT  REAL*8  (A-H,0-Z) 

INTEGER  MCAP,NCAP,INCR,BNDS 
REAL*8MACH,G,D,PHIL,l,SIGMA,Ro,Ao,XSIG,YSIG,PI, 

♦  lw,r2,RAD,ZovL(90),X(90,9 1 ),  Y(90,91  ),Z(90,9 1 ) 

CHARACTER*21  OUTPUT 

COMMON  MCAP,NCAP,PI,INCR,BNDS,PHIL 

C  SET  CONSTANTS 

PI  =  4.dO*DATAN(l.dO) 

RAD  =  PI/180. 

C  CONSTANTS  WHICH  EFFECT  THE  GEOMETRY  OF  THE  WAVERIDER 
MACH  =  10.0 
D  =  5.5  ♦  RAD 
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PHIL  =  50.0  ♦  RAD 
1  =  88.41675 
TYPE  =  1 


C  BASIC  CALCULATIONS  NEEDED  FOR  WAVERIDER  CROSSPL\NE 
C  DEVELOPMENT 

G=  1.4 
MCAP  =  45 
NCAP  =  90 

INCR  =  (2*MCAP)  +  1 
ENDS  =  1 

OUTPUT  =  ’3dsurf 

SIGMA  =  ((G+1.0)/2.0  +1.0/((MACH*D)**2))**0.5 

XSIG  =  SIGMA*COS(PHIL) 

YSIG  =  SIGMA*SIN(PHIL) 

C  VALUE  OF  TYPE  DEFINES  TYPE  OF  PARABOLIC  WAVERIDER. 
DEFAULT 

C  IS  SET  TO  TANGENT  PARABOLIC 

IF  (TYPE.EQ.O)  THEN 
Ro  =  XSIG/2.0 
Ao  =  Ro/(YSIG)**2 
ELSE 

Ro  =  0.75*XSIG 
Ao  =  0.25*XSIG/(YSIG)**2 
END  IF 

Iw  =  1*(1.0  -  Ro/SIGMA) 
rz  =  l*(Ro/SIGMA) 

C  GET  INPUT  DATA 

C  CALL  INPUT(MACH,DJPHIL,l,TYPE,MCAP,NCAP,OUTPUT) 

C  COMPUTE  WAVERIDER  GEOMETRY 

CALL  WAVEBODY(Ro,Ao,14),SIGMA,rz,Rin,Rcb^vL, 

*  X,YJZ) 


119 


C  CALL  OUTDAT  TO  WRITE  TO  A  DATA  FILE 
CALL  OUTDAT(Rin,Rcb,lw,l.MACH,D^vL, 

*  X,Y^) 

END 

C  END  OF  MAIN  PROGRAM 
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SUBROUTINE  WAVEBODY(Ro,A,l,D,SIGMA,rz,Rin,Rcb, 

*  ZovL,X,Y^) 

C  WAVEBODY  DEFINES  THE  X,Y,&  Z  COORDINATES  OF  CROSSPLANES 
C  OF  THE  WAVERIDER  CONFIGURATION 

IMPLICIT  REAL*8  (A-H,0-Z) 

INTEGER  UK,INCR>1,N,BNDS 
REAL*8PHI(90)^HIZ(90)JUn,rz,ts(90)^vL(90), 

*  Rcb,Xin(90,60),Yin(90,60);Zin(90,60), 

*  Xcb(90,60),Ycb(90,60)^cb(90,60),AJlo,l, 

*  SIGMA,test,X(90.91),Y(90,91), 

*  Z(90,91)ADELTA(90),A1  ,PLANE(90),AOJ»ROG 

COMMON  MCAP,NCAP,PI,INCR,BNDS,PHIL 

C  SET  INITIAL  VALUE  OF  PHI  TO  COINCIDE  WITH  THE 
C  NUMBER  OF  CROSSPLANES:  (PHIL/RAD)/NCAP 

AO  =  .001 

PHI(1)  =  0.0*PI/180. 

PROG  =  PHIL-PHI(l) 

CALL  GEOM(PROG,AO,NCAP-14>LANE) 

DO  N  =  2;^CAP 

PHI(N)  =  PHI(N-1)  +  PLANE(NCAP+1-N) 

END  DO 

PHI(NCAP)  =  PHIL 
scale  =  1*D 
DO  10  I  =  1,NCAP 

C  CORRECT  FOR  ERRORS  IN  MACHINE  ZERO 

test  =  ((COS(PHI(I)))**2-4.0*Ro 

*  *A*(SIN(PHI(I)))**2) 

IF(test.LE.0.0)THEN 
test  =  0.0 
END  IF 

rs(I)  =  (ySIGMA)*(2.’»Roy(COS(PHI(I))+(test)**0.5) 
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ZovL(I)  =  (rs(I)-rz)/(l*(1.0-Ro/SIGMA)) 

C  A1  DETERMINES  THE  INITIAL  PACKING  INCREMENT  OF  BODY  POINTS 
C  FROM  THE  LEADING  EDGE  OUTWARD. 


A1  =  .0005 


C  COMPUTE  THE  PHI  INCREMENTS  FOR  BODY  POINT  GENERATION 
STEP  =  PHI® 

CALL  GEOM  (STEP,A1>1CAP^ELTA) 

PHIZ(l)  =  0.0 
DO  15  J  =  2>ICAP+1 

15  PHIZ®  =  PHIZ(J-1)  +  DELTA(MCAP+2-J) 

PHIZ(MCAP+1)  =  PHI(I) 

DO  20  J  =  1>ICAP+1 
K  =  (I-1)*(MCAP+1)+J 
if  (k  .eq.  1)  then 

Rinl  =  ((2.*  Ro)/(COS(PHIZ(J))+((COS(PHIZ(J)))**2 
*  -  4.*Ro*A*(SIN(PHIZ(J)))**2)**0.5)) 

end  if 


Rin  =  ((2.*  Ro)/(COS(PHlZ(J))+((COS(PHIZ(J)))**2 

♦  -  4.*Ro*A*(SIN(PHIZ(J)))**2)**0.5)) 

Rcb  =  (((ZovL(I)*(1.0-Ro/SIGMA)+Ro/SIGMA)**2  + 

♦  (SIGMA**2  -  L0)*Rin**2/SIGMA**2)**0.5) 

C  CALCULATE  THE  X,Y,&  Z  VALUES  FOR  EACH  CROSSPLANE. 
C  Z  IS  A  CONSTANT  FOR  EACH  CROSSPLANE 

Xin(IJ)  =  Rin*COS(PHIZ(J)) 

Yin(I®  =  Rin*SIN(PHIZ(J)) 

Zin(I,J)  =  rs(I) 

Xcb(I,J)  =  Rcb*COS(PHIZ(J)) 

Ycba,J)  =  Rcb*SIN(PHIZ(J)) 

Zcb®  J)  =  rs® 

20  CONTINUE 
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10  CONTINUE 


C  COMBINE  THE  X,  Y,  &  Z  VALUES  OF  THE  FREESTREAM  AND 
C  COMPRESSION  SURFACES  INTO  ONE  ARRAY 

DO  30  I  =  1,  NCAP 
DO  40  J  =  I>ICAP 


Y(I,J)  =  Xin(I,J)*scale  -  Rinl 

1(1,1)  =  Yin(IJ)*scale 

X(I,J)  =  Zin(I,J)  -  (l*Ro/SIGMA) 

40  CONTINUE 
30  CONTINUE 

DO  50  M  =  1,NCAP 
DO  60  N  =  MCAP+1,INCR 

Y(M,N)  =  Xcb(M,(INCR+l)-N)*scale  -  Rinl 
Z(M,N)  =  Ycb(M,(INCR+l)-N)*scale 
X(M^  =  Zcb(M,(INCR+l)-N)  -  (l*Ro/SIGMA) 
60  CONTINUE 

50  CONTINUE 

RETURN 

END 
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4ci|i  *  4> «  *  *  <l>  4>  *  *  <l>  *  <«>  4>  4>  <t>  I*  <(■  <t>  *  *  :)<  <l>  *  >t<  *  *  *  i|>  *  *  *  «  « >l>  *  *  *  *  <» 


SUBROUTINE  OUTDAT(RinJlcb,lw,l,MACH,D^vL,X,Y^) 

C  OUTPUT  SENDS  RESULTS  TO  AN  OUTPUT  FILE 
c  This  subroutine  also  adds  the  spline  fit  required  for  a 
c  viscous  waverider. 
c 

INTEGER  I,J,INCR,BNDS,nim,dim,  ibcbeg,  ibcl,  nose 
REALMS  X(90,91),Y(90,91),Z(90,91),lw,MACaD^n,Rcb,ZovL(90) 
REAL*8  h,xl  ,x2,x3 ,x4,y  1  ,y2,y3,y4,slope,displ 
REALMS  xsp(20),  ysp(20),  tau(2),  c(4,2),  sep,  xsep 
REAL*8  xnose(20),ynose(40)^ose(20Xyspcb(20) 

COMMON  MCAP,NCAP,PI,INCR,BNDS,PHIL 

OPEN(UNIT=12,file=’visc.dat’,form=’formatted’) 
OPEN(UNIT=14,file=’curve.dat’,form=’formatted’) 
c  write(12,*)  ’NETWORK=WAVERIDER’,mcap*2+19,ncap,  ’  NEW’ 

Print*,  ’Enter  distance  (meters)  to  separate  surfaces’ 

Read*,  sep 

xsep  =  0.5  *  sep 
yl=y(l,l) 
c 
c 

c  split  the  two  surfaces  by  1  cm.  and  put  the  waverider  nose  at  Y=0 

c 

c 

Do  10  j=l,ncap 
Do  20  i=mcap+l,mcap*2+l 
y(j,i)=y(j.i>i-sep 
20  continue 

10  continue 

do  15  j=l,ncq) 
do  25  i=l,mcjq)*2+l 
y(j,i)=y(j,i)-yi-xsep 
25  continue 

15  continue 

nose=l 

xsep=xsep/3.0 
100  fonnat(3(2x,E15.8)) 
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444  foiTnat(8(2x,el  5 . 8),2x,i3) 

445  fonnat(5(2x,el5.8),2x,’Xplane  =  ’,i3) 

446  format(lx’first  row  is  xl-x5,  second  is  yl-y5’) 

Print*,’Enter  slope  for  spline  ends’ 
read*,  slope 

Print*,’Enter  number  of  points  for  spline’ 
read*,  dim 

write(12,*)  ’title='’  "’ 
write(12,*)  ’variables=x,y,z’ 
write( 1 2,222)  mcap*  2- 1  +dim*2,ncap+dim- 1 
c  write(12,222)  mcap*2-l+dim*2 
222  format(lx,’zone  t="  ",  i=’,i4,’,  j=’,i4,’,  f=point’) 
c222  format(lx,’2one  t="  ",  i=’,i4,’,  f=point’) 

c 

c 

c  Set  up  to  calculate  the  splines  for  the  blunted  nose 

c 

c 

zl=sep 

z2=z(2,mcap+l)  +  xsep 

z5=0.0 

y2=y(2,l) 

yl=y(l,l) 

y2=y(2,l) 

y4=y(2,mcap*2+ 1 ) 

y3==y(l,mcap*2+l)+(y4-y(l,mc^*2+l))/z2*zl 

xl=x(l,l)+(x(2,l)-x(l,l))/z2*zl 

x2=x(2,l) 

x5=xl-xsep 

y5=(yi+y3)*o.5 

c 

c 

c  Write  out  the  nose  coordinates 

c 

c 

Do  35  i=l,mcap*2+dim*2-l 
write(12,100)  xl,y5,z5 
35  continue 
c 
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Calculate  the  X-Y  plane  splines 


tau(l)=x5 

tau(2)=xl 

c(l,l)=y5 

c(l,2)=yl 

c(2,l)=-1.0*slope 

c(2,2)=0.0 

ncub=2 

ibcl=l 

call  cubspl  (tau,c,ncub,ibcl,ibcl) 

h  =  (xl-x5)/FL0AT(dim) 

Do  30  ii=l,dini 
xnose(ii)  =  x5  +  h*FLOAT(ii) 
xt  =  h*FLOAT(ii) 

ynose(ii)  =  c(l,l)  +  xt*(c(2,l)+xt*(c(3,l)+xt*c(4,l) 

&  /3.0)/2.0) 

continue 

c(U)=y3 

c(2,l)=slope 

c(2,2Hy3-yl)/(x2-xl) 

call  cubspl  (tau,c,ncub,ibcl,ibcl) 

Do  40  ii=l,dim 
xt  =  h*FLOAT(ii) 

ynose(ii+dim)  =  c(l,l)  +  xt*(c(2,l)+xt*(c(3,l)+xt*c(4,l) 
&  /3.0)/2.0) 

continue 


Do  the  X-Z  Plane  spline 


c(l,l)=z5 

c(l^)=zl 

c(2,l)=2.0*slope 


c(2,2)=(z2-zl  )/(x2-x  1 ) 

call  cubspl  (tau,c,ncub,ibcl,ibcl) 

Do  50  ii=l,dim 
xt  =  h*FLOAT(ii) 

znose(ii)  =  c(l,l)  +  xt’''(c(2,l)+xt*(c(3,l)+xt*c(4,l) 
&  /3.0)/2.0) 

continue 


Now,  do  the  linear  interpolation  to  fill  in  the  new 
Y-Z  planes  with  X  coordinates  given  by  xnose(i) 
First,  though,  we  need  to  get  the  spline  fit  for 
the  Y-Z  plane  at  X(nose). 


nose  =  2 

yl=y(nose,mcap+l)  -  sep 

y2=y(nose,mcap) 

y  3=y(nose,mcap+ 1 ) 

y4=y(nose,mcap+2) 

y5=(y3+yl)*.5 

xl=z(nose,mcap+l) 

x2=z(nose,mcap) 

x3=z(nose,mcap+ 1 ) 

x4=z(nose,nicap+2) 

x5=xl+xsep 

tau(l)  =  xl 

tau(2)  =  x5 

c(l,l)  =  yl 

c(l^)  =  y5 

C(2,l)  =  (yl-y2)/(xl-x2) 
c(2,2)  =  slope 
ncub  =  2 
ibcl  =  1 


call  cubspl  (tau,c,ncub,ibcl,ibcl) 
h  =  (x5-xl)/FLOAT(dim) 


xsp(ii)  =  xl  +  h*FLOAT(ii) 
xt  =  h*FLOAT(ii) 

ysp(ii)  =  c(l,l)  +  xt*(c(2,l)+xt*(c(3,l)+xt*c(4,l) 

&  /3.0)/2.0) 

60  continue 

tau(l)  =  x3 
tau(2)  =  x5 
c(l.l)  =  y3 
c(1.2)  =  y5 

c(2,l)  =  (y3-y4)/(x3-x4) 
c(2,2)  =  -1.0*slope 
ncub  =  2 
ibcl  =  1 

call  cubspl  (tau,c,ncub,ibcl,ibcl) 

h  =  (x5-x3)/FLOAT(dim) 

Do  70  ii=l,dim 
xsp(ii)  =  x3  +  h*FLOAT(ii) 
xt  =  h*FLOAT(ii) 

yspcb(ii)  =  c(l,l)  +  xt*(c(2,l)+xt*(c(3,l)+xt*c(4,l) 
&  /3.0)/2.0) 

70  continue 

c 

c 

c  Now  start  the  interpolation 

c 

c 

Do  80  j=2,dim-l 
c  write(12,*)  ’zone’ 
c 
c 

c  The  freestream  surface 

c 

c 

Do  90  i=l^cap 
y  1  =(y5-ynoseij))/(y5-y(nose,  1 )) 

&  *(y(nose,i)-y(nose,l))  +  ynose(j) 

zl  =(znose(j)/(z(nose,mcjq)+ 1  )))*z(nose,i) 
write(12,100)  x!K>se(j),yl,zl 
90  continue 
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c 

c 

c  The  freestream  spline 

c 

c 

Do  110  i=l,dim-l 
yl=(y5-ynose(j))/(y5-y(nose,l)) 

&  ♦(ysp(i)-y(nose,l))  +  ynose(j) 

zl  =(znose(j)/(z(nose,mcap+ 1  )))*  xsp(i) 
write(12,100)  xnose(j),yl,zl 
110  continue 
c 
c 

c  The  compression  spline 

c 

c 

Do  120  i=dim-l,l,-l 

yl=ynose(j+dim)-(ynose(j+dim)-y5)/(y(nose,mcap*2+l)-y5) 
&  *(y(nose,mcap*2+l)-yspcb(i)) 

zl  =(znose(j)/(z(nose,mcap+l  )))*xsp(i) 
write(12,100)  xnose(j),yl,zl 
120  continue 
c 
c 

c  The  compression  surface 

c 

c 

Do  130  i=mc£q5+l,2*mc^+l 

y  l=ynose(j+dim)-(ynose(j+dim)-y  5)/(y(nose,mcap*2+ 1  )-y  5) 
&  *(y(nose,mcap*2+l)-y(nose,i)) 

zl  =(znose(j)/(z(nose,mcap+ 1  )))*  z(nose,i) 
write(12,100)  xnose(j),yl,zl 
130  continue 
80  continue 
c 
c 

c  Now  do  the  main  portion  of  the  waverider 

c 

c 

write(6,*)  ’Nose  starts  at  plane’,nose 
Do  140  mm=nose,ncap 
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c  write(12,*)  ’zone’ 

print*,  ’Loop  number’ ,mm 

yl=y(mm,mcap+l)  -  sep 

y2=y(mm,mcap) 

y  3=y  (nun,mcap+ 1 ) 

y4=y(mm,mcap+2) 

y5=0/3+yl)*.5 

X 1  =z(mm,mcap+ 1 ) 

x2=z(mm,mcap) 

x3=z(mm,mcap+ 1 ) 

x4=z(mm,mcap+2) 

x5=xl+xsep 

tau(l)  =  xl 

tau(2)  =  x5 

c(l,l)  =  yl 

c(l,2)  =  y5 

c(2,l)  =  (yl-y2)/(xl-x2) 
c(2,2)  =  slope 
ncub  =  2 
ibcl  =  1 


call  cubspl  (tau,c,ncub,ibcl,ibcl) 

h  =  (x5-xl)/FLOAT(dim) 

Do  150  ii=l,dim-l 
xsp(ii)  =  xl  +  h*FLOAT(ii) 
xt  =  h*FLOAT(ii) 

ysp(ii)  =  c(l,l)  +  xt*(c(2,l)+xt*(c(3,l)+xt*c(4,l) 
&  /3.0)/2.0) 

ISO  continue 

Do  160  ii=l,mc^ 

write(12,100)  x(mm,ii),y(mm,ii),z(mm,ii) 

160  0>ntinue 

Do  170  ii=l,dim-l 
write(12,100)  x(nim,ii),ysp(ii),xsp(ii) 

170  continue 

tau(l)  =  x3 
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tau(2)  =  x5 
c(l,l)  =  y3 
c(l,2)  =  y5 

c(2,l)  =  (y3-y4)/(x3-x4) 
c(2,2)  =  -1.0*slope 
ncub  =  2 
ibcl  =  1 


call  cubspl  (tau,c,ncub,ibcl,ibcl) 

h  =  (x5-x3)/FLOAT(dim) 

Do  180  ii=l,dim-l 
xsp(ii)  =  x3  +  h*FLOAT(ii) 
xt  =  h*FLOAT(ii) 

ysp(ii)  =  c(l,l)  +  xt*(c(2,l)+xt*(c(3,l)+xt*c(4,l) 
&  /3.0)/2.0) 

180  continue 

Do  190  ii=dim-l,l,-l 
write(12,100)  x(nim,ii),ysp(ii),xsp(ii) 

190  continue 

Do  200  ii=mcap+l,nicap*2+l 
write(  1 2, 1 00)  x(nim,ii),y(mni,ii),z(mm,ii) 

200  Continue 

140  continue 

111  fonnat(lx,3(el 5.8,2x),’Free  Stream  Surface’) 

112  fonnat(lx,3(el5.8,2x),’  Compression  Surface’) 

113  fonnat(5x,’X’,17x,’Y’,17x,’Z’,17x,’Plane  :  ’,14) 
c 

RETURN 

END 
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SUBROUTINE  GEOM  (sum,al,num,astep) 

IMPLICIT  REALMS  (a-h,o-2) 

REAL*  8  astep(num),xsum,fevaljac,dr,rinit, 
*  rdif,fstop,rmax 

rinit  =  l.OdO 
fstop  =  0.0001 OdO 
icount  =  0 
imax  =  10 
rmax  =  1.5d0 
rdif=0.01d0 
xsum  =  sum/al 
1  rinit  =  rinit  +  rdif 
dr  =  rinit 
icount  =  0 

IF  (rinit.gt.miax)  THEN 
WRITE  (*,*)  ’  rinit  exceeded  rmax’ 

GO  TO  99 
END  IF 

10  CONTINUE 
feval  =  l.OdO 
DO  20  i=l,num-l 
20  feval  =  feval  +  dr**i 
feval  =  feval  -  xsum 
IF  (ABS(feval).lt.fstop)  GO  TO  89 
jac  =  1  .OdO 
DO  30  i=l,num-2 
30  jac  =  jac  +  (i+l)*dr**i 
dr  =  dr  -  feval/jac 
icount  =  icount  +  1 
IF  (icount.eq.imax)  THEN 
GOTO  1 
ELSE 
GOTO  10 
END  IF 

89  astep(l)  =  al 
r  =  dr 

DO  40  i=2,num 
40  astep(i)  =  al*r**(i-l) 
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99  RETURN 
END 

SUBROUTINE  CUBSPL  (TAU,  C,  N,  IBCBEG,  IBCEND) 

^  ifi^HfiifiHi*****************  ******************************** 

C  N=Number  of  data  points.  Assumed  to  be  >=2. 
c  (tau(i),  c(l,i),  i=l,n)=Abcissa  and  ordinates  of  the 
c  data  points.  Tau  is  assumed  to  be  strictly  increasing, 

c  IBCBEG,  IBCEND  =  boundary  condition  indicators  and 

c  c(2,l),  c(2,n)  =  boundary  condition  information,  specifically, 
c  IBCBEG=0  means  no  boundary  condition  at  tau)l)  is  given, 
c  In  this  case,  the  not-a-knot  condition  is  used,  i.e.  the 

c  jump  in  the  third  derivative  across  tau(2)  is  forced  to 

c  0,  thus  the  Erst  and  second  cubic  polynomial  pieces 

c  are  made  to  coincide. 

c  IBCBEG=1  means  that  the  slope  at  tau(l)  is  made  to  equal 
c  to  c(2,l),  supplied  by  input. 

c  IBCBEG=2  means  that  the  second  derivative  at  tau(l)  is 
c  made  equal  to  c(2,l),  supplied  by  input, 

c  IBCEND=0,  1,  or  2  has  analogous  meaning  concerning  the 
c  boundary  condition  at  tau(n),  with  the  additional 

information  taken  from  c(2,n). 

OUTPUT  ******************************* 
c(j,i),  j=l,  4;  1=1,  L  (L=N-1)  =  the  polynomial  coefTiecients 
of  the  cubic  interpolating  spline  with  the  interior  knots  (or 
joints)  tau(2),  tau(n-l).  Precisely,  in  the  interval  (tau(I), 
tau(I+l)),  the  spline  F  is  given  by 
F(x)=c(I,l)+h*(c(2,I)+h*(c(3,I)+h*c(4,I)/3)/2) 
where  h=x-tau(I).  The  function  program  *PPVALU*  may  be 
used  to  evaluate  F  or  its  derivatives  from  tau,  c,  L  and  K=4 

INTEGER  IBCBEG,  IBCEND,N,I,J,L,M 
REALMS  C(4,N),TAU(N),DIVDF1,DIVDF343TAU,G 
C 

L  =  N-1 
C 

DO  10  M=2,N 

C(3>1)  =  TAU(M)  -  TAU(M-1) 

C(4>1)  =  (C(1,M)  -  C(1,M-1))/C(3,M) 

10  CONTINUE 
C 

ff(IBCBEG-l)  11,15,16 
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11  IF  (N  .GT.2)  GOTO  12 
C 

C(4,l)=1.0 
C(3,l)  =  1.0 
C(2,l)  =  2.0*C(4,2) 

GOTO  25 
C 

12  C(4,l)  =  C(3,3) 

C(3,l)  =  C(3,2)  +  C(3,3) 

C(2’l)  =  ((C(3,2)+2.0*C(3.1))*C(4.2)*C(3,3)+C(3,2)**2*C(4,3)) 

&  /C(3,l) 

GOTO  19 
C 

15  C(4,l)=1.0 
C(3,l)  =  0.0 
GOTO  18 

C 

16  C(4,l)  =  2.0 
C(3,l)=  1.0 

C(2,l)  =  3.0*C(4,2)-C(3,2)/2.0*C(2,1) 

18  IF  (N  .EQ.  2)  GOTO  25 
C 

19  DO  20 

G  =  -C(3,M+1)/C(4,M-1) 

C(2,M)  =  G*C(2>I-1)+3.0*(C(3,M)*C(4,M+1)+C(3,M+1)*C(4,M+1)) 
C(4,M)  =  G*C(3^-1)+2.0*(C(3,M)+C(3,M+1)) 

20  CONTINUE 
C 

IF  (IBCEND  -1)  21,30,24 

21  IF  (N  .EQ.  3  .AND.  IBCBEG  EQ.  0)  GOTO  22 
C 

G  =  C(3,N-1)+C(3,N) 

C(2,N)  =  ((C(3,N)+2.0*G)*C(4,N)*C(3,N-1) 

&  +C(3,N)**2*(C(l,N-l)-C(l,N-2))/C(3,N-l))/G 

G  =  -G/C(4,N-1) 

C(4,N)  =  C(3,N-1) 

GOTO  29 
C 

22  C(2,N)  =  2.0*C(4,N) 

C(4,N)=  I.O 
GOTO  28 

C 

24  C(2,N)  =  3.0*C(4,N)+C(3,N)/2.0*C(2,N) 
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C(4,N)  =  2.0 
GOTO  28 

25  IF  (IBCEND  -  1)  26,30,24 

26  IF  (IBCBEG  GT.  0)  GOTO  22 
C 

C(2,N)  =  C(4,N) 

GOTO  30 

28  G  =  -1.0/C(4,N-1) 

C 

29  C(44^  =  G*C(3,N-1)  +  C(4,N) 

C(2,N)  =  (G*C(2,N-1)+C(2,N))/C(4,N) 

C 

30  DO  40  J=L,1,-1 

C(2,J)  =  (C(2,J)-C(3,J)*C(2.J+1))/C(4.J) 

40  CONTINUE 
C 

DO  50  I=2,N 
DTAU  =  C(3,I) 

DIVDFl  =  (C(1,I)-C(1,I-1))/DTAU 
DIVDF3  =  C(2,M)+C(2,I)-2.0*DIVDF1 
C(3,I-1)  =  2.0*(DIVDF1-C(2.I-1)-DIVDF3)/DTAU 
C(4,I-1)  =  (DIVDF3/DTAU)*(6.0/DTAU) 

50  CONTINUE 

RETURN 

END 


135 


APPENDIX  D:  CNIDAT  CONTROL  FILE 


ICASE,IMTHDX,IMTEIDY,IMTHDZ  l=Roe,2=McC,3=SW,4=vL 

22  1  1  1 

NEND 
100 
INS 
1 

IL,  JL,  KL 
30,  47,  53 

ILCTST,ICFL,CFLEXP,CFLMAX,CFL 

1  5  1  0.9  0.01 

IREST,CFCRHO,CFCEI,CFLPEN,CEXPPEN,INOFRZ 
1  0.1  0.1  1.  0.  5 

(I/J/K)LMTR,OMEGA,OETA(I/J/K),R(I/J/K)DELT,  (I/J/K)ENTH,(I/J/K)ISO 

2  2  2  1.  -1  -1  -1  0.05  0.05  0.05  2  2  2  0  1  1 

R(I/J/K)SPR 

0  0  0 

IDD(l/2/3/4/5) 

0  10  0  1 
JDD(l/2/3/4/5) 

0  0  10  1 
KDD(l/2/3/4/5) 

000  1  1 

IPC,IMPLT,NSWPS,COEF,AA 
2  12  0.  1.5 

IADBWL,ALPHA4»HI,TWALL,RMINF,REL,  RLSCL,TINT,IMETRC 
0  0.000000  0.  530.00  10.  1.0653E6  1.00  408.57  0 

ITURB,ITBFRZ,APLUS,CCP,  CKLEB,CWK,  SMLK,BIGK,  PRT, 

CMUTHISIMTR 

0  1  26.  2.08  0.3  0.25  0.4  1.68E-02  0.9  14.  0 

IREAD,IGRID,IP3DOP,IDGBUG>IODPR,IP3DMD,NRST,IFMRTI,IFMRTO,IINT 

0  2  2  0  10  500000  0  1  1  1 

lEXRNGS 

0 


IF  lEXRNGS  NON  ZERO  THEN,FOR 

IDMY=1, lEXRNGS 

CHRDMY, 
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IEXMTHD,IIEXRL,IIEXRH,UEXRL,IJEXRH,IKEXRL,IKEXRH,IEXLMT 


JEXRNGS 

0 


IF  JEXRNGS  NON  ZERO  THEN, 

FOR  JDMY=1,JEXRNGS 
CHRDMY 

JEXMTHDJIEXRL,JIEXRH,JJEXRL,JJEXRH,JKEXRL,JKEXRH,JEXLMT 


KEXRNGS 

0 


IF  KEXRNGS  NON  ZERO  THEN,FOR 

KDMY=1,KEXRNGS 

CHRDMY 

KEXMTHD,KIEXRL,KIEXRH,KJEXRL,KJEXRH,KKEXRL,KKEXRH,KEXLMT 


BOUNDARY  CONDITIONS:  Convention 

IBCTYPE=1  =>  Freestream  fixed  boundary  condition 

2  =>  Solid  boundary 

3  — >  Symmetry  boundary 

4  =>  Singular  boundary 

5  =>  Zero  gradient  boundary  conditon 

6  =>  Constant  gradient  boundary  condition 

7  =>  None 

8  =>  Angle  of  attack  (not  mature) 

9  ==>  Pie 

10  =>  Periodic 

1 1  =>  Characteristic  boundary  conditions 
IFACORD(IFACE),IRNGS(IFACE)  Face  1 

1  1 

F  0  RE  A  C  HR  N  G 

IBC(ST/EN),JBC(ST/EN),KBC(ST/EN),IBCTYPE,(W/W)UBC,(I/J/K)INDBC 

1  1  1  47  1  53  1  0  0  0  0  0  0 
Face  2 

2  1 

IBC(ST/EN),JBC(ST/EN),KBC(ST/EN),IBCTYPE,(W/W)UBC,(I/J/K)INDBC 
30  30  1  47  1  53  5  1  1  1  999  999  999 
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Face  3 
3  2 

IBC(ST/EN),JBC(ST/EN),KBC(ST/EN),IBCTYPE,(WAV)UBC.(I/J/K:)INDBC 
1  3  1  1  1  53  9  1  1  1  0  0  1 

IBC(ST/EN),JBC(ST/EN),KBC(ST/EN),IBCTYPE,(WAV)UBC,(I/J/K)INDBC 

3  30  1  1  1  53  2  1  1  1  999  999  999 
Face  4 

4  1 

IBC(ST/EN),JBC(ST/EN),KBC(ST/EN),IBCTYPE.(WAV)UBC,(I/J/K)INDBC 

1  30  47  47  1  53  5  1  1  1  999  999  999 
Face  5 

5  1 

IBC(ST/EN),JBC(ST/EN),KBC(ST/EN),IBCTYPE,(WAV)UBC,(I/J/K)INDBC 

1  30  1  47  1  1  3  1  1  -1  999  999  999 
Face  6 

6  1 

IBC(ST/EN),JBC(ST/EN),KBC(ST/EN),roCTYPE,(WAV)UBC,(I/J/K)INDBC 
1  30  1  47  53  53  3  1  1  -1  999  999  999 


138 


Bibliography 


1.  Anderson,  Dale  A.  et  al.  Computational  Fluid  Mechanics  and  Heat  Transfer. 
New  York,  Hemisphere  Publishing,  1984. 

2.  Anderson,  John  D.  Fundamentals  of  Aerodynamics.  New  York,  McGraw-Hill, 
1984. 

3.  Anderson,  John  D.  Hypersonic  and  High  Temperature  Gas  Dynamics.  New 
York,  McGraw-Hill,  1989. 

4.  Anderson,  John  D.  et  al.  Hypersonic  Waveriders  for  Planetary  Atmospheres. 
AIAA  90-0538  28“'  Aerspace  Sciences  Meeting  and  Exhibit.  Reno  Nevada, 
January  1990. 

5.  Anderson,  John  D.  et  al.  Hypersonic  Waveriders  for  High  Altitude 
Applications.  AIAA  91-0530  29““  Aerspace  Sciences  Meeting  and  Exhibit.  Reno 
Nevada,  January  1991. 

6.  Anderson,  John  D.  et  al.  "Hypersonic  Waveriders:  Effects  of  Chemically 
Reacting  Flow  and  Viscous  Interaction",  AIAA  92-0302  30“’  Aerspace  Sciences 
Meeting  and  Exhibit.  Reno  Nevada,  January  1992. 

7.  Anderson,  John  D.  Modern  Compressible  Flow  With  Historical  Perspective. 
New  York,  McGraw-Hill,  1982. 

8.  Beran,  P.  Class  handout  distributed  in  Aero  752,  Computational  Fluid 
Dynamics.  School  of  Engineering,  Air  Force  Institute  of  Technology  (AU), 
Wright-Patterson  AFB  Ohio,  March  1992. 

9.  Beran  P.  Class  handout  distributed  in  Aero  753,  Computational  Fluid 
Dynamics.  School  of  Engineering,  Air  Force  Institute  of  Technology  (AU), 
Wright-Patterson  AFB  Ohio,  July  1992. 

10.  Bowcutt,  K.G.  et  al.  Viscous  Optimized  Hypersonic  Waveriders,  AIAA  87-0272, 
January  1987. 

11.  Bowcutt,  K.G.  et  al.  "Numerical  Optimization  of  Conical  Flow  Waverider 
Including  Detailed  Viscous  Effects",  Aerodynamics  of  Hypersonic  Lifting 
Vehicles,  AGARD-CPP-428.  March  1987 

12.  Corda,  Stephen,  et  al.  Viscous  Optimized  Waveriders,  AIAA  88-0369  26“’ 
Aerospace  Sciences  Meeting  and  E^ibit.  Reno  Nevada,  January  1988. 


139 


13.  Corda,  Stephen.  Viscous  Optimized  Waveriders  Designed  from  Flows  over 
Cones  and  Minimum  Drag  Bodies.  PhD  Dissertation,  Department  of  Aerospace 
Engineering,  University  of  Maryland,  College  Park  Maryland,  1988. 

14.  de  Boor,  Carl.  A  Practical  Guide  to  Splines.  New  York,  Springer-Verlag, 
1978. 

15.  Gaitonde,  Datta.  Computation  of  Viscous  Shock/Shock  Hypersonic  Interactions 
with  an  Implicit  Flux  Split  Scheme.  Final  Report,  1  September  1989-1 
September  1990.  Wright-Patterson  AFB  Ohio:  Universal  Energy  Systems, 
December  1990,  (WRDC-TR-90-3076. 

16.  Gaitonde,  Datta.  The  Performance  of  Flux-Split  Algorithms  in  High-Speed 
Flows.  AIAA  92-0186  30*  Aerospace  Sciences  Meeting  and  Exhibit.  Reno 
Nevada,  January  1992. 

17.  Gaitonde,  Datta.  Personal  Interviews.  WL/FIMM,  Wright-Patterson  AFB 
Ohio,  1  July  through  30  October  1992. 

18.  Hasen,  G.  Class  handout  distributed  in  Aero  624,  Advanced  Hypersonics. 
School  of  Engineering,  Air  Force  Institute  of  Technology  (AU),  Wright- 
Patterson  AFB  Ohio,  July  1992. 

19.  Hayes,  Wallace  D.  and  Probstein,  Ronald  F.  Hypersonic  Flow  Theory^  New 
York,  Academic  Press,  1959. 

20.  Kuchemann,  Dietrich.  The  Aerodynamic  Design  of  Aircraft.  New  York, 
Pergamon  Press,  1978. 

21.  McLauglhin,  Thomas  A.  Viscous  Optimized  Waveriders  for  Chemiatl 
Equilibrium  Flow,  M.S.  Thesis.  Department  of  Aerospace  Engineering, 
University  of  Maryland,  College  Park  Maryland,  1990. 

22.  Muller,  B.  et  al.  "Simulation  of  Hypersonic  Waverider  Flow".  Proceeding  of 
the  1st  International  Hypersonic  Waverider  Symposium.  University  of  Maryland, 
College  Park  Maryland,  October  1990. 

23.  Murthy,  T.  K.  S.  Computational  Methods  in  Hypersonic  Aerodynamics. 
Computational  Mechanics  Publications,  Kluwer  Aerodynamics  Publishers,  1992. 

24.  Nonweiler,  Terence  R.F.  "The  Waverider  Wing  In  Retrospect  and  Prospect  -  A 
Personalized  View",  Proceeding  of  the  1st  International  Hypersonic  Wxenir 
Symposium.  University  of  Maryland,  College  Park  Maryland,  October  1990. 


140 


25.  Rasmussen,  Maurice  L.  Optimization  of  Waverider  Configurations  Generated 
from  Axisymmetric  Conical  Flows. 

26.  Rasmussen,  Maurice  L.  et  al.  On  Waverider  Shapes  Applied  to  Aero-Space 
Plane  Forebody  Configurations. 

27.  Rasmussen,  Maurice  L.  Class  Notes  distributed  in  Aero  630,  School  of 
Engineering,  Air  Force  Institute  of  Technology,  Wright-Patterson  AFB,  Ohio 

28.  Stecklein,  Capt  Gregory  O.  A  Comparative  Study  of  Numerical  Versus 
Analytical  Waverider  Solutions.  M.S.  thesis,  AFrr/ENY/GAE91D-26.  School 
of  Engineering,  Air  Force  Institute  of  Technology  (AU),  Wright-Patterson  AFB 
Ohio,  December  1991. 

29.  Steger,  Joseph  L.  and  Warming,  R.F.  "Flux  Vector  Splitting  of  the  Inviscid 
Gasdynamic  Equations  with  Application  to  Finite-Difference  Methods",  Journal 
of  Computational  Physics,  40:  263-293  (April  1980) 

30.  Steinbrenner,  John  P.  et  al.  The  GRIDGEN  3D  Multiple  Block  Grid  Generation 
System:  Interim  Report,  1  October  1987-1  October  1990.  Contract  F33615-87- 
C-3003.  Fort  Worth  Texas:  General  Dynamics  Corporation,  April  1991 
(WRDC-TR-90-3022). 

31.  Takashima,  N.  et  al.  "Navier-Stokes  Computation  of  a  Viscous  Optimized 
Waverider",  AIAA  92-0305  30*  Aerspace  Sciences  Meeting  and  Exhibit.  Reno 
Nevada,  January  1992. 

32.  Townend,  L.H.  "Research  and  Design  for  Lifting  Reentry",  Progress  in 
Aerospace  Sciences.  18:1-80,  1979. 

33.  Vanmol,  D.  Aerodynamic  Heating  to  Hypersonic  Waveriders,  M.S.  Thesis, 
Department  of  Aerospace  Engineering,  University  of  Maryland,  College  Park 
Maryland,  1991 

34.  Vincenti,  Walter  D.  et  al.  Introduction  to  Physical  Gas  Dynamics.  Malabar 
Florida,  John  Wiley  &  Sons,  1965. 

35.  Wc-T.iing,  R.F  et  al.  Math  Comp.  29.  1975  1037 

36.  White,  Frank  M.  Viscous  Fluid  Flow  (Second  Edition).  New  York,  MGav- 
Hill,  1991. 

37.  Yee,  H.C.  A  Class  of  High-Resolution  Explicit  and  Implicit  Shock-Capturing 


141 


Methods.  NASA  Technical  Memorandum  101088,  Ames  Research  Center, 
California,  February  1989. 


142 


Vita 


First  Lieutenant  James  A.  Mundy  V  was  bom  on  24  July  1965  in  Roanoke 
Virginia.  He  received  his  GED  from  the  state  of  Massachusetts  in  1986,  and  enlisted 
in  the  United  States  Air  Force.  He  attended  the  University  of  Maryland,  graduating  with 
a  Bachelor  of  Science  degree  in  Aerospace  Engineering  in  1988.  Upon  graduation  he 
attended  Air  Force  Officer  Candidate  School  and  received  a  reserve  commission  in  the 
USAF.  He  began  as  a  test  engineer  for  the  49*  Test  and  Evaluation  Squadron  at 
Barksdale  AFB,  Louisianna,  where  he  was  involved  in  the  operational  teseting  and 
evaluation  of  Strategic  Air  Command’s  airborne  weapons  systems.  Lt  Mundy  served 
at  the  49*  Test  Squadron  until  May  of  1991,  when  he  entered  the  School  of 
Engineering,  Air  Force  Institute  of  Technology. 


143 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
0MB  No.  0704-0188 


Puolic  reporttng  sauroen  ror  this  vCUectJOO  of  information  ts  estimated  to  average  t  nour  oer  response,  including  the  time  tor  reviewing  mstruaicns.  searching  OKisting  oata  sources, 
gathering  jnd  maintaining  the  oata  needed,  and  completing  and  reviewing  the  colleaion  of  information.  Send  comments  reaarding  this  burden  estimate  or  any  other  aspea  of  this 
coMection  of  information,  including  suggestions  for  reducing  this  Durden,  to  Washington  HeadQuarters  Services.  Directorate  tor  information  Operations  and  Reports.  121S  Jefferson 
Davis  Highway.  Suite  1204.  Arlington.  VA  22202*4302.  and  to  the  Office  of  Management  and  Budget.  Paperwork  Reduction  Projea  (0704-0188).  Washington.  DC  20503. 


1.  AGENCY  USE  ONLY  (Leave  blank) 


4.  TITLE  AND  SUBTITLE 


2.  REPORT  DATE 

Decenber  1992 


3.  REPORT  TYPE  AND  OATES  COVERED 

Master's  Thesis 


The  Effects  of  Viscosity  on  a  Conically 
Derived  Haverider 


6.  AUTHOR(S) 

Jaaies  A.  Mundy,  ILt,  USAF 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  AOORESS(E$) 

Air  Force  Institute  of  Technology 
Hright-Patterson  AFB,  Ohio  45433 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

AFIT/6AE/EMT/92D-23 


9.  SPONSORING /MONITORING  AGENCY  NAME(S)  AND  AODRESS{ES) 

WL/FI 

Wright-Patterson  AFB,  Ohio  45433 


10.  SPONSORING /MONITORING 
AGENCY  REPORT  NUMBER 


12a.  DISTRIBUTION /AVAILABILITY  STATEMENT 


12b.  DISTRIBUTION  CODE 


Unlinited  Distribution 


13.  ABSTRACT  (Maximum  200  words) 

This  study  investigated  the  effects  of  the  interaction  between  the  viscous 
boundary  layer  and  the  shock  wave  produced  by  a  Mach  10  inviscid  optimized 
waverlder.  An  iiQillcit,  Roe  flux-splitting  algorithm,  developed  by  MZ./FIMN,  was 
used  to  solve  the  flow  field.  A  validation  for  the  inviscid  version  of  the  CFD 
algorithm  was  accoiqilished  by  cosparing  the  numerical  data  produced  by  the  CFD 
code  to  the  analytic  results  derived  by  Rasmussen,  and  by  coaparlson  to  results 
of  the  explicit  version  of  the  same  Roe  flux-splitting  code.  The  coaputational 
results  compared  favorably.  The  Inviscid  case  studied  using  the  ixqplicit  code 
produced  results  identical,  for  all  practical  purposes,  to  those  of  the  explicit 
code,  though  approximately  twice  as  quickly.  The  results  of  the  viscous  flow 
case  matched  well  with  the  results  predicted  by  theory.  The  lift  to  drag  ratio 
calculated,  5.74,  Is  coi^arable  to  the  results  of  other  researchers. 


14.  SUBJEa  TERMS 

Inviscid  Optimized  Waverider,  Navier-Stokes  Solution, 
Viscous  Flow,  Hypersonic  Vehicle 


15.  NUMBER  OF  PAGES 
158 


16.  PRICE  CODE 


17.  SECURITY  CLASSIFICATION 
OF  REPORT 

URCLA8SXFXBD 


18.  SECURITY  CLASSIFICATION 
OF  THIS  PAGE 

UHCXJLSSXFXED 


19.  SECURITY  CLASSIFICATION 
OF  ABSTRACT 

DMCLASSIFIBD 


20.  UMITATION  OF  ABSTRACT 


NSN  7S40-01-280-S500 


Standard  Form  298  (Rev.  2-89) 
Pmcnbcd  by  ansi  Std  ZJS-tS 
2M.ia2 


